From 031f9c018b6a6f1c92e09b010c512e7fe22fbf83 Mon Sep 17 00:00:00 2001 From: Yifan Li Date: Wed, 10 Sep 2025 22:49:57 -0400 Subject: [PATCH 1/3] add fix cboamd for cavity Born-Oppenheimer MD --- source/lmp/fix_cboamd.cpp | 636 +++++++++++++++++++++++++++++ source/lmp/fix_cboamd.h | 147 +++++++ source/lmp/plugin/deepmdplugin.cpp | 12 + 3 files changed, 795 insertions(+) create mode 100644 source/lmp/fix_cboamd.cpp create mode 100644 source/lmp/fix_cboamd.h diff --git a/source/lmp/fix_cboamd.cpp b/source/lmp/fix_cboamd.cpp new file mode 100644 index 0000000000..6df6d04988 --- /dev/null +++ b/source/lmp/fix_cboamd.cpp @@ -0,0 +1,636 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_cboamd.h" + +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "neighbor.h" +#include "update.h" +#include "output.h" +#include "respa.h" +#include "universe.h" +#include "utils.h" + +#include +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixCBOAMD::FixCBOAMD(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), + model_potential(nullptr), + model_dipole(nullptr), + model_polar(nullptr), + deepmd_dipole(nullptr), + deepmd_polar(nullptr), + type_map(nullptr), + ntypes(0), + dipole(nullptr), + polarizability(nullptr), + forces_deepmd(nullptr), + photons_enabled(false), + nphoton(0), + omega_photon(nullptr), + lambda_photon(nullptr), + lambda_vector(nullptr), + qa(nullptr), + pa(nullptr), + fa(nullptr), + ea(nullptr), + dt(0.0), + current_step(0), + output_file(nullptr) +{ + vector_flag = 1; + extvector = 1; + size_vector = 3; + + if (narg < 4) error->all(FLERR,"Illegal fix cboamd command"); + + // Parse command line arguments + int iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg], "dipole") == 0) { + if (iarg + 1 >= narg) error->all(FLERR,"Illegal fix cboamd command"); + model_dipole = utils::strdup(arg[iarg + 1]); + iarg += 2; + } else if (strcmp(arg[iarg], "polar") == 0) { + if (iarg + 1 >= narg) error->all(FLERR,"Illegal fix cboamd command"); + model_polar = utils::strdup(arg[iarg + 1]); + iarg += 2; + } else if (strcmp(arg[iarg], "photons") == 0) { + if (iarg + 1 >= narg) error->all(FLERR,"Illegal fix cboamd command"); + photons_enabled = utils::logical(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "nphoton") == 0) { + if (iarg + 1 >= narg) error->all(FLERR,"Illegal fix cboamd command"); + nphoton = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "omega") == 0) { + if (iarg + 1 + nphoton > narg) error->all(FLERR,"Illegal fix cboamd command"); + memory->create(omega_photon, nphoton, "fix_cboamd:omega_photon"); + for (int i = 0; i < nphoton; i++) { + omega_photon[i] = utils::numeric(FLERR, arg[iarg + 1 + i], false, lmp); + } + iarg += 1 + nphoton; + } else if (strcmp(arg[iarg], "lambda") == 0) { + if (iarg + 1 + nphoton > narg) error->all(FLERR,"Illegal fix cboamd command"); + memory->create(lambda_photon, nphoton, "fix_cboamd:lambda_photon"); + for (int i = 0; i < nphoton; i++) { + lambda_photon[i] = utils::numeric(FLERR, arg[iarg + 1 + i], false, lmp); + } + iarg += 1 + nphoton; + } else if (strcmp(arg[iarg], "lambda_vector") == 0) { + if (iarg + 1 + 3 * nphoton > narg) error->all(FLERR,"Illegal fix cboamd command"); + memory->create(lambda_vector, nphoton, 3, "fix_cboamd:lambda_vector"); + for (int i = 0; i < nphoton; i++) { + for (int j = 0; j < 3; j++) { + lambda_vector[i][j] = utils::numeric(FLERR, arg[iarg + 1 + 3 * i + j], false, lmp); + } + } + iarg += 1 + 3 * nphoton; + } else { + error->all(FLERR,"Unknown fix cboamd keyword: {}", arg[iarg]); + } + } + + // Check required parameters + if (!model_dipole) error->all(FLERR,"fix cboamd requires dipole model"); + // if (!model_polar) error->all(FLERR,"fix cboamd requires polarizability model"); + if (photons_enabled && nphoton <= 0) error->all(FLERR,"fix cboamd: nphoton must be > 0 when photons enabled"); + if (photons_enabled && !omega_photon) error->all(FLERR,"fix cboamd: omega must be specified when photons enabled"); + if (photons_enabled && !lambda_photon) error->all(FLERR,"fix cboamd: lambda must be specified when photons enabled"); + if (photons_enabled && !lambda_vector) error->all(FLERR,"fix cboamd: lambda_vector must be specified when photons enabled"); + // if (dt <= 0.0) error->all(FLERR,"fix cboamd: dt must be > 0"); + + // Set up arrays + ntypes = atom->ntypes; + memory->create(type_map, ntypes + 1, "fix_cboamd:type_map"); + for (int i = 1; i <= ntypes; i++) { + type_map[i] = i - 1; // Default mapping: type 1 -> 0, type 2 -> 1, etc. + } + + // Allocate arrays + memory->create(dipole, 3, "fix_cboamd:dipole"); + memory->create(polarizability, 9, "fix_cboamd:polarizability"); + memory->create(forces_deepmd, atom->nmax * 3, "fix_cboamd:forces_deepmd"); + + if (photons_enabled) { + memory->create(qa, nphoton, "fix_cboamd:qa"); + memory->create(pa, nphoton, "fix_cboamd:pa"); + memory->create(fa, nphoton, "fix_cboamd:fa"); + memory->create(ea, nphoton, "fix_cboamd:ea"); + + // Initialize photon coordinates and momenta + for (int i = 0; i < nphoton; i++) { + qa[i] = 0.0; + pa[i] = 0.0; + fa[i] = 0.0; + ea[i] = 0.0; + } + } + + // Initialize DeepMD models + init_deepmd_models(); + + // Open output file + output_file = fopen("cboamd_output.dat", "w"); + if (!output_file) error->all(FLERR,"Cannot open cboamd output file"); + fprintf(output_file, "# Step Time Energy Dipole_x Dipole_y Dipole_z Pol_xx Pol_xy Pol_xz Pol_yx Pol_yy Pol_yz Pol_zx Pol_zy Pol_zz"); + if (photons_enabled) { + for (int i = 0; i < nphoton; i++) { + fprintf(output_file, " qa_%d pa_%d ea_%d fa_%d", i, i, i, i); + } + } + fprintf(output_file, "\n"); +} + +/* ---------------------------------------------------------------------- */ + +FixCBOAMD::~FixCBOAMD() +{ + cleanup_deepmd_models(); + + if (model_potential) delete[] model_potential; + if (model_dipole) delete[] model_dipole; + if (model_polar) delete[] model_polar; + + memory->destroy(type_map); + memory->destroy(dipole); + memory->destroy(polarizability); + memory->destroy(forces_deepmd); + + if (photons_enabled) { + memory->destroy(omega_photon); + memory->destroy(lambda_photon); + memory->destroy(lambda_vector); + memory->destroy(qa); + memory->destroy(pa); + memory->destroy(fa); + memory->destroy(ea); + } + + if (output_file) fclose(output_file); +} + +/* ---------------------------------------------------------------------- */ + +int FixCBOAMD::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= POST_FORCE; + mask |= POST_FORCE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::init() +{ + dt = update->dt; + // Check if we have the right atom style + if (!atom->tag_enable) + error->all(FLERR,"fix cboamd requires atom IDs"); + + // Set up neighbor list + neighbor->add_request(this); +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::setup(int vflag) +{ + post_force(vflag); + // // Initialize DeepMD calculations + // convert_coordinates_to_deepmd_format(); + // compute_deepmd_dipole(); + // // compute_deepmd_polarizability(); + + // if (photons_enabled) { + // compute_cboa_forces(); + // } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::initial_integrate(int /*vflag*/) +{ + for(int alpha = 0; alpha < nphoton; alpha++) { + // Update photon coordinates and momenta using Velocity Verlet + pa[alpha] += 0.5 * fa[alpha] * dt * PS_TO_AU; + qa[alpha] += pa[alpha] * dt * PS_TO_AU; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::final_integrate() +{ + for(int alpha = 0; alpha < nphoton; alpha++) { + // Update photon momenta using Velocity Verlet + pa[alpha] += 0.5 * fa[alpha] * dt * PS_TO_AU; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::post_force(int vflag) +{ + // Update photon coordinates if enabled + // if (photons_enabled) { + // update_photon_coordinates(); + // } + + // Compute DeepMD properties + convert_coordinates_to_deepmd_format(); + compute_deepmd_dipole(); + // compute_deepmd_polarizability(); + + // Compute CBOA forces if photons enabled + if (photons_enabled) { + compute_cboa_forces(); + } + + // Write output + write_output(); + + current_step++; +} + +/* ---------------------------------------------------------------------- */ + +double FixCBOAMD::compute_scalar() +{ + // Return total energy + return 0.0; // This would need to be implemented with actual energy value +} + +/* ---------------------------------------------------------------------- */ + +double FixCBOAMD::compute_vector(int n) +{ + // Return dipole components + // if (n >= 0 && n < 3) return dipole[n]; + for (int i=0; i= 0 && i < 3 && j >= 0 && j < 3) { + return polarizability[i * 3 + j]; + } + return 0.0; +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::write_restart(FILE *fp) +{ + // Write restart data + if (comm->me == 0) { + int size = 0; + if (photons_enabled) size += 4 * nphoton; + fwrite(&size, sizeof(int), 1, fp); + + if (photons_enabled) { + fwrite(qa, sizeof(double), nphoton, fp); + fwrite(pa, sizeof(double), nphoton, fp); + fwrite(fa, sizeof(double), nphoton, fp); + fwrite(ea, sizeof(double), nphoton, fp); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::restart(char *buf) +{ + // Read restart data + int size; + memcpy(&size, buf, sizeof(int)); + buf += sizeof(int); + + if (photons_enabled && size >= 4 * nphoton) { + memcpy(qa, buf, nphoton * sizeof(double)); + buf += nphoton * sizeof(double); + memcpy(pa, buf, nphoton * sizeof(double)); + buf += nphoton * sizeof(double); + memcpy(fa, buf, nphoton * sizeof(double)); + buf += nphoton * sizeof(double); + memcpy(ea, buf, nphoton * sizeof(double)); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::grow_arrays(int nmax) +{ + memory->grow(forces_deepmd, nmax * 3, "fix_cboamd:forces_deepmd"); +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::copy_arrays(int i, int j, int delflag) +{ + // Copy forces for atom i to atom j + for (int k = 0; k < 3; k++) { + forces_deepmd[j * 3 + k] = forces_deepmd[i * 3 + k]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::set_arrays(int i) +{ + // Set forces for atom i + for (int k = 0; k < 3; k++) { + forces_deepmd[i * 3 + k] = 0.0; + } +} + +/* ---------------------------------------------------------------------- */ + +int FixCBOAMD::pack_exchange(int i, double *buf) +{ + // Pack exchange data + int n = 0; + for (int k = 0; k < 3; k++) { + buf[n++] = forces_deepmd[i * 3 + k]; + } + return n; +} + +/* ---------------------------------------------------------------------- */ + +int FixCBOAMD::unpack_exchange(int nlocal, double *buf) +{ + // Unpack exchange data + int n = 0; + for (int k = 0; k < 3; k++) { + forces_deepmd[nlocal * 3 + k] = buf[n++]; + } + return n; +} + +/* ---------------------------------------------------------------------- */ + +int FixCBOAMD::pack_restart(int i, double *buf) +{ + // Pack restart data + int n = 0; + for (int k = 0; k < 3; k++) { + buf[n++] = forces_deepmd[i * 3 + k]; + } + return n; +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::unpack_restart(int nlocal, int nth) +{ + // Unpack restart data + // This method is called when restarting from a checkpoint + // For now, we'll just initialize photon coordinates to zero + if (photons_enabled) { + for (int i = 0; i < nphoton; i++) { + qa[i] = 0.0; + pa[i] = 0.0; + fa[i] = 0.0; + ea[i] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixCBOAMD::size_restart(int nlocal) +{ + return 3; +} + +/* ---------------------------------------------------------------------- */ + +int FixCBOAMD::maxsize_restart() +{ + return 3; +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::init_deepmd_models() +{ + try { + // Initialize DeepMD models using C++ interface + deepmd_dipole = new deepmd_compat::DeepTensor(model_dipole); + if (model_polar) deepmd_polar = new deepmd_compat::DeepTensor(model_polar); + + if (comm->me == 0) { + utils::logmesg(lmp, "DeepMD models initialized successfully:\n"); + // utils::logmesg(lmp, " Potential: {}\n", model_potential); + utils::logmesg(lmp, " Dipole: {}\n", model_dipole); + if (model_polar) utils::logmesg(lmp, " Polarizability: {}\n", model_polar); + } + } catch (const std::exception& e) { + error->all(FLERR, "Failed to initialize DeepMD models: {}", e.what()); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::cleanup_deepmd_models() +{ + // Clean up DeepMD models + if (deepmd_dipole) { + delete deepmd_dipole; + deepmd_dipole = nullptr; + } + if (deepmd_polar) { + delete deepmd_polar; + deepmd_polar = nullptr; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::compute_deepmd_dipole() +{ + try { + // Compute dipole using DeepMD + + // dipole_grad_deepmd is the negative gradient of dipole w.r.t. atomic coordinates + // It is a flattened array of size 3 (dipole components) * N (atoms) * 3 (atomic coordinate components) + deepmd_dipole->compute(dipole_deepmd, dipole_grad_deepmd, + dipole_virial_deepmd, dipole_atom_deepmd, dipole_atom_virial_deepmd, + coords_deepmd, atom_types_deepmd, cell_deepmd); + // Extract dipole components (DeepMD returns in eV/A, convert to a.u.) + // dipole[0] = dipole_deepmd[0] * ANGSTROM_TO_BOHR; + // dipole[1] = dipole_deepmd[1] * ANGSTROM_TO_BOHR; + // dipole[2] = dipole_deepmd[2] * ANGSTROM_TO_BOHR; + dipole[0] = dipole_deepmd[0]; + dipole[1] = dipole_deepmd[1]; + dipole[2] = dipole_deepmd[2]; + } catch (const std::exception& e) { + error->all(FLERR, "DeepMD dipole computation failed: {}", e.what()); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::compute_deepmd_polarizability() +{ + try { + // Compute polarizability using DeepMD + deepmd_polar->compute(polar_deepmd, polar_grad_deepmd, polar_virial_deepmd, + polar_atom_deepmd, polar_atom_virial_deepmd, coords_deepmd, atom_types_deepmd, cell_deepmd); + + // Extract polarizability components (DeepMD returns in eV/A, convert to a.u.) + for (int i = 0; i < 9; i++) { + polarizability[i] = polar_deepmd[i] * EV_PER_ANGSTROM_TO_HARTREE_PER_BOHR; + } + } catch (const std::exception& e) { + error->all(FLERR, "DeepMD polarizability computation failed: {}", e.what()); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::update_photon_coordinates() +{ + // Update photon coordinates using velocity Verlet + for (int i = 0; i < nphoton; i++) { + // Update position + qa[i] += (pa[i] + fa[i] * dt * 0.5) * dt; + + // Store old force + double fa_old = fa[i]; + + // Update force (this would be computed based on dipole and polarizability) + fa[i] = -ea[i] * omega_photon[i]; + + // Update momentum + pa[i] += (fa_old + fa[i]) * dt * 0.5; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::compute_cboa_forces() +{ + // Compute cavity Born-Oppenheimer forces + // This would implement the CBOA force calculation + // For now, we'll just set dummy values + + for (int i = 0; i < nphoton; i++) { + ea[i] = omega_photon[i] * qa[i]; + for (int j = 0; j < 3; j++) { + ea[i] -= lambda_photon[i] * lambda_vector[i][j] * dipole[j]; + } + } + + int nlocal = atom->nlocal; + double **f = atom->f; + + for (int dp = 0; dp < 3; dp++) { // dp is the direction of dipole or photon + for (int i = 0; i < nlocal; i++) { + for (int di = 0; di < 3; di++) { // di is the direction of atomic coordinate + int idx = dp * nlocal * 3 + i * 3 + di; // index in dipole_grad_deepmd: dipole component dp, atom i, coordinate di + for (int alpha = 0; alpha < nphoton; alpha++) { + // CBOA force contribution from photon alpha + f[i][di] -= ea[alpha] * lambda_photon[alpha] * lambda_vector[alpha][di] * dipole_grad_deepmd[idx] / EV_TO_HARTREE; + } + } + } + } + + for (int alpha = 0; alpha < nphoton; alpha++) { + fa[alpha] = -ea[alpha] * omega_photon[alpha]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::write_output() +{ + // Write output to file + if (comm->me == 0) { + fprintf(output_file, "%d %.6f %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e", + current_step, update->dt * current_step, + dipole[0], dipole[1], dipole[2], + polarizability[0], polarizability[1], polarizability[2], + polarizability[3], polarizability[4], polarizability[5], + polarizability[6], polarizability[7], polarizability[8]); + + if (photons_enabled) { + for (int i = 0; i < nphoton; i++) { + fprintf(output_file, " %.6e %.6e %.6e %.6e", qa[i], pa[i], ea[i], fa[i]); + } + } + fprintf(output_file, "\n"); + fflush(output_file); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCBOAMD::convert_coordinates_to_deepmd_format() +{ + int nlocal = atom->nlocal; + int *type = atom->type; + double **x = atom->x; + double *boxlo = domain->boxlo; + double *boxhi = domain->boxhi; + + // Resize vectors + coords_deepmd.resize(nlocal * 3); + atom_types_deepmd.resize(nlocal); + cell_deepmd.resize(9); + + // Convert coordinates from LAMMPS to DeepMD format + for (int i = 0; i < nlocal; i++) { + coords_deepmd[i * 3] = x[i][0]; // x + coords_deepmd[i * 3 + 1] = x[i][1]; // y + coords_deepmd[i * 3 + 2] = x[i][2]; // z + atom_types_deepmd[i] = type_map[type[i]]; + } + + // Set cell (assuming orthorhombic box) + double lx = boxhi[0] - boxlo[0]; + double ly = boxhi[1] - boxlo[1]; + double lz = boxhi[2] - boxlo[2]; + + cell_deepmd[0] = lx; cell_deepmd[1] = 0.0; cell_deepmd[2] = 0.0; + cell_deepmd[3] = 0.0; cell_deepmd[4] = ly; cell_deepmd[5] = 0.0; + cell_deepmd[6] = 0.0; cell_deepmd[7] = 0.0; cell_deepmd[8] = lz; +} diff --git a/source/lmp/fix_cboamd.h b/source/lmp/fix_cboamd.h new file mode 100644 index 0000000000..3d9bc5d413 --- /dev/null +++ b/source/lmp/fix_cboamd.h @@ -0,0 +1,147 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(cboamd,FixCBOAMD); +// clang-format on +#else + +#ifndef LMP_FIX_CBOAMD_H +#define LMP_FIX_CBOAMD_H + +#include "fix.h" + +// DeePMD interface includes +#ifdef DP_USE_CXX_API +#ifdef LMPPLUGIN +#include "DeepPot.h" +#include "DeepTensor.h" +#else +#include "deepmd/DeepPot.h" +#include "deepmd/DeepTensor.h" +#endif +namespace deepmd_compat = deepmd; +#else +#ifdef LMPPLUGIN +#include "deepmd.hpp" +#else +#include "deepmd/deepmd.hpp" +#endif +namespace deepmd_compat = deepmd::hpp; +#endif + +namespace LAMMPS_NS { + +class FixCBOAMD : public Fix { + public: + FixCBOAMD(class LAMMPS *, int, char **); + ~FixCBOAMD() override; + int setmask() override; + void init() override; + void setup(int) override; + void initial_integrate(int) override; + void final_integrate() override; + void post_force(int) override; + double compute_scalar() override; + double compute_vector(int) override; + double compute_array(int, int) override; + void write_restart(FILE *) override; + void restart(char *) override; + void grow_arrays(int) override; + void copy_arrays(int, int, int) override; + void set_arrays(int) override; + int pack_exchange(int, double *) override; + int unpack_exchange(int, double *) override; + int pack_restart(int, double *) override; + void unpack_restart(int, int) override; + int size_restart(int) override; + int maxsize_restart() override; + + private: + // DeepMD model paths + char *model_potential; + char *model_dipole; + char *model_polar; + + // DeepMD model objects using C++ interface + deepmd_compat::DeepTensor *deepmd_dipole; + deepmd_compat::DeepTensor *deepmd_polar; + + // Type mapping for atoms + int *type_map; + int ntypes; + + // Output arrays + double *dipole; + double *polarizability; + double *forces_deepmd; + + // DeepMD computation arrays + std::vector coords_deepmd; + std::vector atom_types_deepmd; + std::vector cell_deepmd; + std::vector dipole_deepmd; + std::vector dipole_grad_deepmd; + std::vector dipole_virial_deepmd; + std::vector dipole_atom_deepmd; + std::vector dipole_atom_virial_deepmd; + std::vector polar_deepmd; + std::vector polar_grad_deepmd; + std::vector polar_virial_deepmd; + std::vector polar_atom_deepmd; + std::vector polar_atom_virial_deepmd; + + // Cavity Born-Oppenheimer parameters + bool photons_enabled; + int nphoton; + double *omega_photon; + double *lambda_photon; + double **lambda_vector; + + // Photon coordinates and momenta + double *qa, *pa, *fa, *ea; + + // Time step + double dt; + + // Current step + int current_step; + + // Output file + FILE *output_file; + + // Helper functions + void init_deepmd_models(); + void cleanup_deepmd_models(); + void compute_deepmd_dipole(); + void compute_deepmd_polarizability(); + void update_photon_coordinates(); + void compute_cboa_forces(); + void write_output(); + + // DeepMD computation helpers + void convert_coordinates_to_deepmd_format(); + + // Unit conversion constants + static constexpr double EV_TO_HARTREE = 0.0367493; + static constexpr double ANGSTROM_TO_BOHR = 1.88973; + static constexpr double EV_PER_ANGSTROM_TO_HARTREE_PER_BOHR = 0.0194467; + static constexpr double PS_TO_AU = 41.3413745758e3; // 1 fs = 41.3413745758 a.u. of time + static constexpr double HARTREE_PER_BOHR_TO_EV_PER_ANGSTROM = 51.42208619083232; // 1 Hartree/Bohr = 51.42208619083232 eV/Angstrom +}; + +} + +#endif +#endif diff --git a/source/lmp/plugin/deepmdplugin.cpp b/source/lmp/plugin/deepmdplugin.cpp index d3b54f8e41..06bbaf68b0 100644 --- a/source/lmp/plugin/deepmdplugin.cpp +++ b/source/lmp/plugin/deepmdplugin.cpp @@ -4,6 +4,7 @@ */ #include "compute_deeptensor_atom.h" #include "deepmd_version.h" +#include "fix_cboamd.h" #include "fix_dplr.h" #include "lammpsplugin.h" #include "pair_deepmd.h" @@ -26,6 +27,10 @@ static Fix* fixdplr(LAMMPS* lmp, int narg, char** arg) { return new FixDPLR(lmp, narg, arg); } +static Fix *fixcboamd(LAMMPS *lmp, int narg, char **arg) { + return new FixCBOAMD(lmp, narg, arg); +} + #if LAMMPS_VERSION_NUMBER >= 20220328 static KSpace* pppmdplr(LAMMPS* lmp) { return new PPPMDPLR(lmp); } #endif @@ -66,6 +71,13 @@ extern "C" void lammpsplugin_init(void* lmp, void* handle, void* regfunc) { plugin.creator.v2 = (lammpsplugin_factory2*)&fixdplr; (*register_plugin)(&plugin, lmp); + plugin.style = "fix"; + plugin.name = "cboamd"; + plugin.info = "fix cboamd " STR_GIT_SUMM; + plugin.author = "DeePMD-kit"; + plugin.creator.v2 = (lammpsplugin_factory2 *)&fixcboamd; + (*register_plugin)(&plugin, lmp); + #if LAMMPS_VERSION_NUMBER >= 20220328 // lammps/lammps# plugin.style = "kspace"; From 2f88d065ab92c6bad5af2a16f04aeb8cfb6e1c94 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 11 Sep 2025 03:21:42 +0000 Subject: [PATCH 2/3] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/lmp/fix_cboamd.cpp | 411 ++++++++++++++++------------- source/lmp/fix_cboamd.h | 69 ++--- source/lmp/plugin/deepmdplugin.cpp | 4 +- 3 files changed, 264 insertions(+), 220 deletions(-) diff --git a/source/lmp/fix_cboamd.cpp b/source/lmp/fix_cboamd.cpp index 6df6d04988..010f94fb9d 100644 --- a/source/lmp/fix_cboamd.cpp +++ b/source/lmp/fix_cboamd.cpp @@ -1,3 +1,4 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,6 +14,10 @@ #include "fix_cboamd.h" +#include +#include +#include + #include "atom.h" #include "atom_vec.h" #include "comm.h" @@ -22,126 +27,153 @@ #include "memory.h" #include "modify.h" #include "neighbor.h" -#include "update.h" #include "output.h" #include "respa.h" #include "universe.h" +#include "update.h" #include "utils.h" -#include -#include -#include - using namespace LAMMPS_NS; using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixCBOAMD::FixCBOAMD(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - model_potential(nullptr), - model_dipole(nullptr), - model_polar(nullptr), - deepmd_dipole(nullptr), - deepmd_polar(nullptr), - type_map(nullptr), - ntypes(0), - dipole(nullptr), - polarizability(nullptr), - forces_deepmd(nullptr), - photons_enabled(false), - nphoton(0), - omega_photon(nullptr), - lambda_photon(nullptr), - lambda_vector(nullptr), - qa(nullptr), - pa(nullptr), - fa(nullptr), - ea(nullptr), - dt(0.0), - current_step(0), - output_file(nullptr) -{ +FixCBOAMD::FixCBOAMD(LAMMPS* lmp, int narg, char** arg) + : Fix(lmp, narg, arg), + model_potential(nullptr), + model_dipole(nullptr), + model_polar(nullptr), + deepmd_dipole(nullptr), + deepmd_polar(nullptr), + type_map(nullptr), + ntypes(0), + dipole(nullptr), + polarizability(nullptr), + forces_deepmd(nullptr), + photons_enabled(false), + nphoton(0), + omega_photon(nullptr), + lambda_photon(nullptr), + lambda_vector(nullptr), + qa(nullptr), + pa(nullptr), + fa(nullptr), + ea(nullptr), + dt(0.0), + current_step(0), + output_file(nullptr) { vector_flag = 1; extvector = 1; size_vector = 3; - if (narg < 4) error->all(FLERR,"Illegal fix cboamd command"); - + if (narg < 4) { + error->all(FLERR, "Illegal fix cboamd command"); + } + // Parse command line arguments int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg], "dipole") == 0) { - if (iarg + 1 >= narg) error->all(FLERR,"Illegal fix cboamd command"); + if (iarg + 1 >= narg) { + error->all(FLERR, "Illegal fix cboamd command"); + } model_dipole = utils::strdup(arg[iarg + 1]); iarg += 2; } else if (strcmp(arg[iarg], "polar") == 0) { - if (iarg + 1 >= narg) error->all(FLERR,"Illegal fix cboamd command"); + if (iarg + 1 >= narg) { + error->all(FLERR, "Illegal fix cboamd command"); + } model_polar = utils::strdup(arg[iarg + 1]); iarg += 2; } else if (strcmp(arg[iarg], "photons") == 0) { - if (iarg + 1 >= narg) error->all(FLERR,"Illegal fix cboamd command"); + if (iarg + 1 >= narg) { + error->all(FLERR, "Illegal fix cboamd command"); + } photons_enabled = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; } else if (strcmp(arg[iarg], "nphoton") == 0) { - if (iarg + 1 >= narg) error->all(FLERR,"Illegal fix cboamd command"); + if (iarg + 1 >= narg) { + error->all(FLERR, "Illegal fix cboamd command"); + } nphoton = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; } else if (strcmp(arg[iarg], "omega") == 0) { - if (iarg + 1 + nphoton > narg) error->all(FLERR,"Illegal fix cboamd command"); + if (iarg + 1 + nphoton > narg) { + error->all(FLERR, "Illegal fix cboamd command"); + } memory->create(omega_photon, nphoton, "fix_cboamd:omega_photon"); for (int i = 0; i < nphoton; i++) { omega_photon[i] = utils::numeric(FLERR, arg[iarg + 1 + i], false, lmp); } iarg += 1 + nphoton; } else if (strcmp(arg[iarg], "lambda") == 0) { - if (iarg + 1 + nphoton > narg) error->all(FLERR,"Illegal fix cboamd command"); + if (iarg + 1 + nphoton > narg) { + error->all(FLERR, "Illegal fix cboamd command"); + } memory->create(lambda_photon, nphoton, "fix_cboamd:lambda_photon"); for (int i = 0; i < nphoton; i++) { lambda_photon[i] = utils::numeric(FLERR, arg[iarg + 1 + i], false, lmp); } iarg += 1 + nphoton; } else if (strcmp(arg[iarg], "lambda_vector") == 0) { - if (iarg + 1 + 3 * nphoton > narg) error->all(FLERR,"Illegal fix cboamd command"); + if (iarg + 1 + 3 * nphoton > narg) { + error->all(FLERR, "Illegal fix cboamd command"); + } memory->create(lambda_vector, nphoton, 3, "fix_cboamd:lambda_vector"); for (int i = 0; i < nphoton; i++) { for (int j = 0; j < 3; j++) { - lambda_vector[i][j] = utils::numeric(FLERR, arg[iarg + 1 + 3 * i + j], false, lmp); + lambda_vector[i][j] = + utils::numeric(FLERR, arg[iarg + 1 + 3 * i + j], false, lmp); } } iarg += 1 + 3 * nphoton; } else { - error->all(FLERR,"Unknown fix cboamd keyword: {}", arg[iarg]); + error->all(FLERR, "Unknown fix cboamd keyword: {}", arg[iarg]); } } - + // Check required parameters - if (!model_dipole) error->all(FLERR,"fix cboamd requires dipole model"); - // if (!model_polar) error->all(FLERR,"fix cboamd requires polarizability model"); - if (photons_enabled && nphoton <= 0) error->all(FLERR,"fix cboamd: nphoton must be > 0 when photons enabled"); - if (photons_enabled && !omega_photon) error->all(FLERR,"fix cboamd: omega must be specified when photons enabled"); - if (photons_enabled && !lambda_photon) error->all(FLERR,"fix cboamd: lambda must be specified when photons enabled"); - if (photons_enabled && !lambda_vector) error->all(FLERR,"fix cboamd: lambda_vector must be specified when photons enabled"); + if (!model_dipole) { + error->all(FLERR, "fix cboamd requires dipole model"); + } + // if (!model_polar) error->all(FLERR,"fix cboamd requires polarizability + // model"); + if (photons_enabled && nphoton <= 0) { + error->all(FLERR, "fix cboamd: nphoton must be > 0 when photons enabled"); + } + if (photons_enabled && !omega_photon) { + error->all(FLERR, + "fix cboamd: omega must be specified when photons enabled"); + } + if (photons_enabled && !lambda_photon) { + error->all(FLERR, + "fix cboamd: lambda must be specified when photons enabled"); + } + if (photons_enabled && !lambda_vector) { + error->all( + FLERR, + "fix cboamd: lambda_vector must be specified when photons enabled"); + } // if (dt <= 0.0) error->all(FLERR,"fix cboamd: dt must be > 0"); - + // Set up arrays ntypes = atom->ntypes; memory->create(type_map, ntypes + 1, "fix_cboamd:type_map"); for (int i = 1; i <= ntypes; i++) { type_map[i] = i - 1; // Default mapping: type 1 -> 0, type 2 -> 1, etc. } - + // Allocate arrays memory->create(dipole, 3, "fix_cboamd:dipole"); memory->create(polarizability, 9, "fix_cboamd:polarizability"); memory->create(forces_deepmd, atom->nmax * 3, "fix_cboamd:forces_deepmd"); - + if (photons_enabled) { memory->create(qa, nphoton, "fix_cboamd:qa"); memory->create(pa, nphoton, "fix_cboamd:pa"); memory->create(fa, nphoton, "fix_cboamd:fa"); memory->create(ea, nphoton, "fix_cboamd:ea"); - + // Initialize photon coordinates and momenta for (int i = 0; i < nphoton; i++) { qa[i] = 0.0; @@ -150,14 +182,18 @@ FixCBOAMD::FixCBOAMD(LAMMPS *lmp, int narg, char **arg) : ea[i] = 0.0; } } - + // Initialize DeepMD models init_deepmd_models(); - + // Open output file output_file = fopen("cboamd_output.dat", "w"); - if (!output_file) error->all(FLERR,"Cannot open cboamd output file"); - fprintf(output_file, "# Step Time Energy Dipole_x Dipole_y Dipole_z Pol_xx Pol_xy Pol_xz Pol_yx Pol_yy Pol_yz Pol_zx Pol_zy Pol_zz"); + if (!output_file) { + error->all(FLERR, "Cannot open cboamd output file"); + } + fprintf(output_file, + "# Step Time Energy Dipole_x Dipole_y Dipole_z Pol_xx Pol_xy Pol_xz " + "Pol_yx Pol_yy Pol_yz Pol_zx Pol_zy Pol_zz"); if (photons_enabled) { for (int i = 0; i < nphoton; i++) { fprintf(output_file, " qa_%d pa_%d ea_%d fa_%d", i, i, i, i); @@ -168,19 +204,24 @@ FixCBOAMD::FixCBOAMD(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -FixCBOAMD::~FixCBOAMD() -{ +FixCBOAMD::~FixCBOAMD() { cleanup_deepmd_models(); - - if (model_potential) delete[] model_potential; - if (model_dipole) delete[] model_dipole; - if (model_polar) delete[] model_polar; - + + if (model_potential) { + delete[] model_potential; + } + if (model_dipole) { + delete[] model_dipole; + } + if (model_polar) { + delete[] model_polar; + } + memory->destroy(type_map); memory->destroy(dipole); memory->destroy(polarizability); memory->destroy(forces_deepmd); - + if (photons_enabled) { memory->destroy(omega_photon); memory->destroy(lambda_photon); @@ -190,14 +231,15 @@ FixCBOAMD::~FixCBOAMD() memory->destroy(fa); memory->destroy(ea); } - - if (output_file) fclose(output_file); + + if (output_file) { + fclose(output_file); + } } /* ---------------------------------------------------------------------- */ -int FixCBOAMD::setmask() -{ +int FixCBOAMD::setmask() { int mask = 0; mask |= INITIAL_INTEGRATE; mask |= FINAL_INTEGRATE; @@ -208,27 +250,26 @@ int FixCBOAMD::setmask() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::init() -{ +void FixCBOAMD::init() { dt = update->dt; // Check if we have the right atom style - if (!atom->tag_enable) - error->all(FLERR,"fix cboamd requires atom IDs"); - + if (!atom->tag_enable) { + error->all(FLERR, "fix cboamd requires atom IDs"); + } + // Set up neighbor list neighbor->add_request(this); } /* ---------------------------------------------------------------------- */ -void FixCBOAMD::setup(int vflag) -{ +void FixCBOAMD::setup(int vflag) { post_force(vflag); // // Initialize DeepMD calculations // convert_coordinates_to_deepmd_format(); // compute_deepmd_dipole(); // // compute_deepmd_polarizability(); - + // if (photons_enabled) { // compute_cboa_forces(); // } @@ -236,9 +277,8 @@ void FixCBOAMD::setup(int vflag) /* ---------------------------------------------------------------------- */ -void FixCBOAMD::initial_integrate(int /*vflag*/) -{ - for(int alpha = 0; alpha < nphoton; alpha++) { +void FixCBOAMD::initial_integrate(int /*vflag*/) { + for (int alpha = 0; alpha < nphoton; alpha++) { // Update photon coordinates and momenta using Velocity Verlet pa[alpha] += 0.5 * fa[alpha] * dt * PS_TO_AU; qa[alpha] += pa[alpha] * dt * PS_TO_AU; @@ -247,9 +287,8 @@ void FixCBOAMD::initial_integrate(int /*vflag*/) /* ---------------------------------------------------------------------- */ -void FixCBOAMD::final_integrate() -{ - for(int alpha = 0; alpha < nphoton; alpha++) { +void FixCBOAMD::final_integrate() { + for (int alpha = 0; alpha < nphoton; alpha++) { // Update photon momenta using Velocity Verlet pa[alpha] += 0.5 * fa[alpha] * dt * PS_TO_AU; } @@ -257,18 +296,17 @@ void FixCBOAMD::final_integrate() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::post_force(int vflag) -{ +void FixCBOAMD::post_force(int vflag) { // Update photon coordinates if enabled // if (photons_enabled) { - // update_photon_coordinates(); + // update_photon_coordinates(); // } - + // Compute DeepMD properties convert_coordinates_to_deepmd_format(); compute_deepmd_dipole(); // compute_deepmd_polarizability(); - + // Compute CBOA forces if photons enabled if (photons_enabled) { compute_cboa_forces(); @@ -276,28 +314,32 @@ void FixCBOAMD::post_force(int vflag) // Write output write_output(); - + current_step++; } /* ---------------------------------------------------------------------- */ -double FixCBOAMD::compute_scalar() -{ +double FixCBOAMD::compute_scalar() { // Return total energy - return 0.0; // This would need to be implemented with actual energy value + return 0.0; // This would need to be implemented with actual energy value } /* ---------------------------------------------------------------------- */ -double FixCBOAMD::compute_vector(int n) -{ +double FixCBOAMD::compute_vector(int n) { // Return dipole components // if (n >= 0 && n < 3) return dipole[n]; - for (int i=0; i= 0 && i < 3 && j >= 0 && j < 3) { return polarizability[i * 3 + j]; @@ -316,14 +357,15 @@ double FixCBOAMD::compute_array(int i, int j) /* ---------------------------------------------------------------------- */ -void FixCBOAMD::write_restart(FILE *fp) -{ +void FixCBOAMD::write_restart(FILE* fp) { // Write restart data if (comm->me == 0) { int size = 0; - if (photons_enabled) size += 4 * nphoton; + if (photons_enabled) { + size += 4 * nphoton; + } fwrite(&size, sizeof(int), 1, fp); - + if (photons_enabled) { fwrite(qa, sizeof(double), nphoton, fp); fwrite(pa, sizeof(double), nphoton, fp); @@ -335,13 +377,12 @@ void FixCBOAMD::write_restart(FILE *fp) /* ---------------------------------------------------------------------- */ -void FixCBOAMD::restart(char *buf) -{ +void FixCBOAMD::restart(char* buf) { // Read restart data int size; memcpy(&size, buf, sizeof(int)); buf += sizeof(int); - + if (photons_enabled && size >= 4 * nphoton) { memcpy(qa, buf, nphoton * sizeof(double)); buf += nphoton * sizeof(double); @@ -355,15 +396,13 @@ void FixCBOAMD::restart(char *buf) /* ---------------------------------------------------------------------- */ -void FixCBOAMD::grow_arrays(int nmax) -{ +void FixCBOAMD::grow_arrays(int nmax) { memory->grow(forces_deepmd, nmax * 3, "fix_cboamd:forces_deepmd"); } /* ---------------------------------------------------------------------- */ -void FixCBOAMD::copy_arrays(int i, int j, int delflag) -{ +void FixCBOAMD::copy_arrays(int i, int j, int delflag) { // Copy forces for atom i to atom j for (int k = 0; k < 3; k++) { forces_deepmd[j * 3 + k] = forces_deepmd[i * 3 + k]; @@ -372,8 +411,7 @@ void FixCBOAMD::copy_arrays(int i, int j, int delflag) /* ---------------------------------------------------------------------- */ -void FixCBOAMD::set_arrays(int i) -{ +void FixCBOAMD::set_arrays(int i) { // Set forces for atom i for (int k = 0; k < 3; k++) { forces_deepmd[i * 3 + k] = 0.0; @@ -382,8 +420,7 @@ void FixCBOAMD::set_arrays(int i) /* ---------------------------------------------------------------------- */ -int FixCBOAMD::pack_exchange(int i, double *buf) -{ +int FixCBOAMD::pack_exchange(int i, double* buf) { // Pack exchange data int n = 0; for (int k = 0; k < 3; k++) { @@ -394,8 +431,7 @@ int FixCBOAMD::pack_exchange(int i, double *buf) /* ---------------------------------------------------------------------- */ -int FixCBOAMD::unpack_exchange(int nlocal, double *buf) -{ +int FixCBOAMD::unpack_exchange(int nlocal, double* buf) { // Unpack exchange data int n = 0; for (int k = 0; k < 3; k++) { @@ -406,8 +442,7 @@ int FixCBOAMD::unpack_exchange(int nlocal, double *buf) /* ---------------------------------------------------------------------- */ -int FixCBOAMD::pack_restart(int i, double *buf) -{ +int FixCBOAMD::pack_restart(int i, double* buf) { // Pack restart data int n = 0; for (int k = 0; k < 3; k++) { @@ -418,8 +453,7 @@ int FixCBOAMD::pack_restart(int i, double *buf) /* ---------------------------------------------------------------------- */ -void FixCBOAMD::unpack_restart(int nlocal, int nth) -{ +void FixCBOAMD::unpack_restart(int nlocal, int nth) { // Unpack restart data // This method is called when restarting from a checkpoint // For now, we'll just initialize photon coordinates to zero @@ -435,32 +469,29 @@ void FixCBOAMD::unpack_restart(int nlocal, int nth) /* ---------------------------------------------------------------------- */ -int FixCBOAMD::size_restart(int nlocal) -{ - return 3; -} +int FixCBOAMD::size_restart(int nlocal) { return 3; } /* ---------------------------------------------------------------------- */ -int FixCBOAMD::maxsize_restart() -{ - return 3; -} +int FixCBOAMD::maxsize_restart() { return 3; } /* ---------------------------------------------------------------------- */ -void FixCBOAMD::init_deepmd_models() -{ +void FixCBOAMD::init_deepmd_models() { try { // Initialize DeepMD models using C++ interface deepmd_dipole = new deepmd_compat::DeepTensor(model_dipole); - if (model_polar) deepmd_polar = new deepmd_compat::DeepTensor(model_polar); - + if (model_polar) { + deepmd_polar = new deepmd_compat::DeepTensor(model_polar); + } + if (comm->me == 0) { utils::logmesg(lmp, "DeepMD models initialized successfully:\n"); // utils::logmesg(lmp, " Potential: {}\n", model_potential); utils::logmesg(lmp, " Dipole: {}\n", model_dipole); - if (model_polar) utils::logmesg(lmp, " Polarizability: {}\n", model_polar); + if (model_polar) { + utils::logmesg(lmp, " Polarizability: {}\n", model_polar); + } } } catch (const std::exception& e) { error->all(FLERR, "Failed to initialize DeepMD models: {}", e.what()); @@ -469,8 +500,7 @@ void FixCBOAMD::init_deepmd_models() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::cleanup_deepmd_models() -{ +void FixCBOAMD::cleanup_deepmd_models() { // Clean up DeepMD models if (deepmd_dipole) { delete deepmd_dipole; @@ -484,16 +514,17 @@ void FixCBOAMD::cleanup_deepmd_models() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::compute_deepmd_dipole() -{ +void FixCBOAMD::compute_deepmd_dipole() { try { // Compute dipole using DeepMD - // dipole_grad_deepmd is the negative gradient of dipole w.r.t. atomic coordinates - // It is a flattened array of size 3 (dipole components) * N (atoms) * 3 (atomic coordinate components) - deepmd_dipole->compute(dipole_deepmd, dipole_grad_deepmd, - dipole_virial_deepmd, dipole_atom_deepmd, dipole_atom_virial_deepmd, - coords_deepmd, atom_types_deepmd, cell_deepmd); + // dipole_grad_deepmd is the negative gradient of dipole w.r.t. atomic + // coordinates It is a flattened array of size 3 (dipole components) * N + // (atoms) * 3 (atomic coordinate components) + deepmd_dipole->compute(dipole_deepmd, dipole_grad_deepmd, + dipole_virial_deepmd, dipole_atom_deepmd, + dipole_atom_virial_deepmd, coords_deepmd, + atom_types_deepmd, cell_deepmd); // Extract dipole components (DeepMD returns in eV/A, convert to a.u.) // dipole[0] = dipole_deepmd[0] * ANGSTROM_TO_BOHR; // dipole[1] = dipole_deepmd[1] * ANGSTROM_TO_BOHR; @@ -508,14 +539,15 @@ void FixCBOAMD::compute_deepmd_dipole() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::compute_deepmd_polarizability() -{ +void FixCBOAMD::compute_deepmd_polarizability() { try { // Compute polarizability using DeepMD - deepmd_polar->compute(polar_deepmd, polar_grad_deepmd, polar_virial_deepmd, - polar_atom_deepmd, polar_atom_virial_deepmd, coords_deepmd, atom_types_deepmd, cell_deepmd); - - // Extract polarizability components (DeepMD returns in eV/A, convert to a.u.) + deepmd_polar->compute(polar_deepmd, polar_grad_deepmd, polar_virial_deepmd, + polar_atom_deepmd, polar_atom_virial_deepmd, + coords_deepmd, atom_types_deepmd, cell_deepmd); + + // Extract polarizability components (DeepMD returns in eV/A, convert to + // a.u.) for (int i = 0; i < 9; i++) { polarizability[i] = polar_deepmd[i] * EV_PER_ANGSTROM_TO_HARTREE_PER_BOHR; } @@ -526,19 +558,18 @@ void FixCBOAMD::compute_deepmd_polarizability() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::update_photon_coordinates() -{ +void FixCBOAMD::update_photon_coordinates() { // Update photon coordinates using velocity Verlet for (int i = 0; i < nphoton; i++) { // Update position qa[i] += (pa[i] + fa[i] * dt * 0.5) * dt; - + // Store old force double fa_old = fa[i]; - + // Update force (this would be computed based on dipole and polarizability) fa[i] = -ea[i] * omega_photon[i]; - + // Update momentum pa[i] += (fa_old + fa[i]) * dt * 0.5; } @@ -546,8 +577,7 @@ void FixCBOAMD::update_photon_coordinates() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::compute_cboa_forces() -{ +void FixCBOAMD::compute_cboa_forces() { // Compute cavity Born-Oppenheimer forces // This would implement the CBOA force calculation // For now, we'll just set dummy values @@ -560,15 +590,20 @@ void FixCBOAMD::compute_cboa_forces() } int nlocal = atom->nlocal; - double **f = atom->f; - - for (int dp = 0; dp < 3; dp++) { // dp is the direction of dipole or photon + double** f = atom->f; + + for (int dp = 0; dp < 3; dp++) { // dp is the direction of dipole or photon for (int i = 0; i < nlocal; i++) { - for (int di = 0; di < 3; di++) { // di is the direction of atomic coordinate - int idx = dp * nlocal * 3 + i * 3 + di; // index in dipole_grad_deepmd: dipole component dp, atom i, coordinate di + for (int di = 0; di < 3; + di++) { // di is the direction of atomic coordinate + int idx = dp * nlocal * 3 + i * 3 + + di; // index in dipole_grad_deepmd: dipole component dp, atom + // i, coordinate di for (int alpha = 0; alpha < nphoton; alpha++) { // CBOA force contribution from photon alpha - f[i][di] -= ea[alpha] * lambda_photon[alpha] * lambda_vector[alpha][di] * dipole_grad_deepmd[idx] / EV_TO_HARTREE; + f[i][di] -= ea[alpha] * lambda_photon[alpha] * + lambda_vector[alpha][di] * dipole_grad_deepmd[idx] / + EV_TO_HARTREE; } } } @@ -581,20 +616,21 @@ void FixCBOAMD::compute_cboa_forces() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::write_output() -{ +void FixCBOAMD::write_output() { // Write output to file if (comm->me == 0) { - fprintf(output_file, "%d %.6f %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e", - current_step, update->dt * current_step, - dipole[0], dipole[1], dipole[2], - polarizability[0], polarizability[1], polarizability[2], + fprintf(output_file, + "%d %.6f %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e " + "%.6e %.6e", + current_step, update->dt * current_step, dipole[0], dipole[1], + dipole[2], polarizability[0], polarizability[1], polarizability[2], polarizability[3], polarizability[4], polarizability[5], polarizability[6], polarizability[7], polarizability[8]); - + if (photons_enabled) { for (int i = 0; i < nphoton; i++) { - fprintf(output_file, " %.6e %.6e %.6e %.6e", qa[i], pa[i], ea[i], fa[i]); + fprintf(output_file, " %.6e %.6e %.6e %.6e", qa[i], pa[i], ea[i], + fa[i]); } } fprintf(output_file, "\n"); @@ -604,33 +640,38 @@ void FixCBOAMD::write_output() /* ---------------------------------------------------------------------- */ -void FixCBOAMD::convert_coordinates_to_deepmd_format() -{ +void FixCBOAMD::convert_coordinates_to_deepmd_format() { int nlocal = atom->nlocal; - int *type = atom->type; - double **x = atom->x; - double *boxlo = domain->boxlo; - double *boxhi = domain->boxhi; - + int* type = atom->type; + double** x = atom->x; + double* boxlo = domain->boxlo; + double* boxhi = domain->boxhi; + // Resize vectors coords_deepmd.resize(nlocal * 3); atom_types_deepmd.resize(nlocal); cell_deepmd.resize(9); - + // Convert coordinates from LAMMPS to DeepMD format for (int i = 0; i < nlocal; i++) { - coords_deepmd[i * 3] = x[i][0]; // x - coords_deepmd[i * 3 + 1] = x[i][1]; // y - coords_deepmd[i * 3 + 2] = x[i][2]; // z + coords_deepmd[i * 3] = x[i][0]; // x + coords_deepmd[i * 3 + 1] = x[i][1]; // y + coords_deepmd[i * 3 + 2] = x[i][2]; // z atom_types_deepmd[i] = type_map[type[i]]; } - + // Set cell (assuming orthorhombic box) double lx = boxhi[0] - boxlo[0]; double ly = boxhi[1] - boxlo[1]; double lz = boxhi[2] - boxlo[2]; - - cell_deepmd[0] = lx; cell_deepmd[1] = 0.0; cell_deepmd[2] = 0.0; - cell_deepmd[3] = 0.0; cell_deepmd[4] = ly; cell_deepmd[5] = 0.0; - cell_deepmd[6] = 0.0; cell_deepmd[7] = 0.0; cell_deepmd[8] = lz; + + cell_deepmd[0] = lx; + cell_deepmd[1] = 0.0; + cell_deepmd[2] = 0.0; + cell_deepmd[3] = 0.0; + cell_deepmd[4] = ly; + cell_deepmd[5] = 0.0; + cell_deepmd[6] = 0.0; + cell_deepmd[7] = 0.0; + cell_deepmd[8] = lz; } diff --git a/source/lmp/fix_cboamd.h b/source/lmp/fix_cboamd.h index 3d9bc5d413..eaa997267f 100644 --- a/source/lmp/fix_cboamd.h +++ b/source/lmp/fix_cboamd.h @@ -1,3 +1,4 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -45,7 +46,7 @@ namespace LAMMPS_NS { class FixCBOAMD : public Fix { public: - FixCBOAMD(class LAMMPS *, int, char **); + FixCBOAMD(class LAMMPS*, int, char**); ~FixCBOAMD() override; int setmask() override; void init() override; @@ -56,37 +57,37 @@ class FixCBOAMD : public Fix { double compute_scalar() override; double compute_vector(int) override; double compute_array(int, int) override; - void write_restart(FILE *) override; - void restart(char *) override; + void write_restart(FILE*) override; + void restart(char*) override; void grow_arrays(int) override; void copy_arrays(int, int, int) override; void set_arrays(int) override; - int pack_exchange(int, double *) override; - int unpack_exchange(int, double *) override; - int pack_restart(int, double *) override; + int pack_exchange(int, double*) override; + int unpack_exchange(int, double*) override; + int pack_restart(int, double*) override; void unpack_restart(int, int) override; int size_restart(int) override; int maxsize_restart() override; private: // DeepMD model paths - char *model_potential; - char *model_dipole; - char *model_polar; - + char* model_potential; + char* model_dipole; + char* model_polar; + // DeepMD model objects using C++ interface - deepmd_compat::DeepTensor *deepmd_dipole; - deepmd_compat::DeepTensor *deepmd_polar; - + deepmd_compat::DeepTensor* deepmd_dipole; + deepmd_compat::DeepTensor* deepmd_polar; + // Type mapping for atoms - int *type_map; + int* type_map; int ntypes; - + // Output arrays - double *dipole; - double *polarizability; - double *forces_deepmd; - + double* dipole; + double* polarizability; + double* forces_deepmd; + // DeepMD computation arrays std::vector coords_deepmd; std::vector atom_types_deepmd; @@ -105,22 +106,22 @@ class FixCBOAMD : public Fix { // Cavity Born-Oppenheimer parameters bool photons_enabled; int nphoton; - double *omega_photon; - double *lambda_photon; - double **lambda_vector; - + double* omega_photon; + double* lambda_photon; + double** lambda_vector; + // Photon coordinates and momenta double *qa, *pa, *fa, *ea; - + // Time step double dt; - + // Current step int current_step; - + // Output file - FILE *output_file; - + FILE* output_file; + // Helper functions void init_deepmd_models(); void cleanup_deepmd_models(); @@ -129,19 +130,21 @@ class FixCBOAMD : public Fix { void update_photon_coordinates(); void compute_cboa_forces(); void write_output(); - + // DeepMD computation helpers void convert_coordinates_to_deepmd_format(); - + // Unit conversion constants static constexpr double EV_TO_HARTREE = 0.0367493; static constexpr double ANGSTROM_TO_BOHR = 1.88973; static constexpr double EV_PER_ANGSTROM_TO_HARTREE_PER_BOHR = 0.0194467; - static constexpr double PS_TO_AU = 41.3413745758e3; // 1 fs = 41.3413745758 a.u. of time - static constexpr double HARTREE_PER_BOHR_TO_EV_PER_ANGSTROM = 51.42208619083232; // 1 Hartree/Bohr = 51.42208619083232 eV/Angstrom + static constexpr double PS_TO_AU = + 41.3413745758e3; // 1 fs = 41.3413745758 a.u. of time + static constexpr double HARTREE_PER_BOHR_TO_EV_PER_ANGSTROM = + 51.42208619083232; // 1 Hartree/Bohr = 51.42208619083232 eV/Angstrom }; -} +} // namespace LAMMPS_NS #endif #endif diff --git a/source/lmp/plugin/deepmdplugin.cpp b/source/lmp/plugin/deepmdplugin.cpp index 06bbaf68b0..bc2be8e955 100644 --- a/source/lmp/plugin/deepmdplugin.cpp +++ b/source/lmp/plugin/deepmdplugin.cpp @@ -27,7 +27,7 @@ static Fix* fixdplr(LAMMPS* lmp, int narg, char** arg) { return new FixDPLR(lmp, narg, arg); } -static Fix *fixcboamd(LAMMPS *lmp, int narg, char **arg) { +static Fix* fixcboamd(LAMMPS* lmp, int narg, char** arg) { return new FixCBOAMD(lmp, narg, arg); } @@ -75,7 +75,7 @@ extern "C" void lammpsplugin_init(void* lmp, void* handle, void* regfunc) { plugin.name = "cboamd"; plugin.info = "fix cboamd " STR_GIT_SUMM; plugin.author = "DeePMD-kit"; - plugin.creator.v2 = (lammpsplugin_factory2 *)&fixcboamd; + plugin.creator.v2 = (lammpsplugin_factory2*)&fixcboamd; (*register_plugin)(&plugin, lmp); #if LAMMPS_VERSION_NUMBER >= 20220328 From 61c7289793b6dd58ab7028f833832f9cdfda7a17 Mon Sep 17 00:00:00 2001 From: Yifan Li Date: Thu, 11 Sep 2025 11:40:57 -0400 Subject: [PATCH 3/3] delete potential-model-related lines --- source/lmp/fix_cboamd.cpp | 5 ----- source/lmp/fix_cboamd.h | 1 - 2 files changed, 6 deletions(-) diff --git a/source/lmp/fix_cboamd.cpp b/source/lmp/fix_cboamd.cpp index 010f94fb9d..c18d3fed4e 100644 --- a/source/lmp/fix_cboamd.cpp +++ b/source/lmp/fix_cboamd.cpp @@ -40,7 +40,6 @@ using namespace FixConst; FixCBOAMD::FixCBOAMD(LAMMPS* lmp, int narg, char** arg) : Fix(lmp, narg, arg), - model_potential(nullptr), model_dipole(nullptr), model_polar(nullptr), deepmd_dipole(nullptr), @@ -207,9 +206,6 @@ FixCBOAMD::FixCBOAMD(LAMMPS* lmp, int narg, char** arg) FixCBOAMD::~FixCBOAMD() { cleanup_deepmd_models(); - if (model_potential) { - delete[] model_potential; - } if (model_dipole) { delete[] model_dipole; } @@ -487,7 +483,6 @@ void FixCBOAMD::init_deepmd_models() { if (comm->me == 0) { utils::logmesg(lmp, "DeepMD models initialized successfully:\n"); - // utils::logmesg(lmp, " Potential: {}\n", model_potential); utils::logmesg(lmp, " Dipole: {}\n", model_dipole); if (model_polar) { utils::logmesg(lmp, " Polarizability: {}\n", model_polar); diff --git a/source/lmp/fix_cboamd.h b/source/lmp/fix_cboamd.h index eaa997267f..bdec51bc0e 100644 --- a/source/lmp/fix_cboamd.h +++ b/source/lmp/fix_cboamd.h @@ -71,7 +71,6 @@ class FixCBOAMD : public Fix { private: // DeepMD model paths - char* model_potential; char* model_dipole; char* model_polar;