From 235372d6e811a229d45d1913b77e413ec3d58827 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Fri, 28 Apr 2023 14:40:59 -0600 Subject: [PATCH 01/41] Fuse some Kokkos kernels to reduce launch latency for small systems --- src/KOKKOS/fix_nve_kokkos.cpp | 30 ++++++++++++++++++++++++ src/KOKKOS/fix_nve_kokkos.h | 20 ++++++++++++++++ src/KOKKOS/modify_kokkos.cpp | 19 +++++++++++++++ src/KOKKOS/modify_kokkos.h | 1 + src/KOKKOS/pair_kokkos.h | 26 +++++++++++++++++++++ src/KOKKOS/verlet_kokkos.cpp | 44 ++++++++++++++++++++++++++++++++--- src/KOKKOS/verlet_kokkos.h | 3 +++ src/fix.cpp | 6 ++--- src/fix.h | 4 +++- src/modify.cpp | 30 ++++++++++++++++++++++++ src/modify.h | 3 +++ src/pair.cpp | 7 +++--- src/pair.h | 3 ++- src/timer.h | 1 + 14 files changed, 185 insertions(+), 12 deletions(-) diff --git a/src/KOKKOS/fix_nve_kokkos.cpp b/src/KOKKOS/fix_nve_kokkos.cpp index 5dcb611d41..f812c575bc 100644 --- a/src/KOKKOS/fix_nve_kokkos.cpp +++ b/src/KOKKOS/fix_nve_kokkos.cpp @@ -29,6 +29,7 @@ FixNVEKokkos::FixNVEKokkos(LAMMPS *lmp, int narg, char **arg) : FixNVE(lmp, narg, arg) { kokkosable = 1; + fuse_integrate_flag = 1; atomKK = (AtomKokkos *) atom; execution_space = ExecutionSpaceFromDevice::space; @@ -168,6 +169,35 @@ void FixNVEKokkos::cleanup_copy() vatom = nullptr; } +/* ---------------------------------------------------------------------- + allow for both per-type and per-atom mass +------------------------------------------------------------------------- */ + +template +void FixNVEKokkos::fused_integrate() +{ + atomKK->sync(execution_space,datamask_read); + atomKK->modified(execution_space,datamask_modify); + + x = atomKK->k_x.view(); + v = atomKK->k_v.view(); + f = atomKK->k_f.view(); + rmass = atomKK->k_rmass.view(); + mass = atomKK->k_mass.view(); + type = atomKK->k_type.view(); + mask = atomKK->k_mask.view(); + int nlocal = atomKK->nlocal; + if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; + + if (rmass.data()) { + FixNVEKokkosFusedIntegrateFunctor functor(this); + Kokkos::parallel_for(nlocal,functor); + } else { + FixNVEKokkosFusedIntegrateFunctor functor(this); + Kokkos::parallel_for(nlocal,functor); + } +} + namespace LAMMPS_NS { template class FixNVEKokkos; #ifdef LMP_KOKKOS_GPU diff --git a/src/KOKKOS/fix_nve_kokkos.h b/src/KOKKOS/fix_nve_kokkos.h index 215aacb4a0..a10f846cba 100644 --- a/src/KOKKOS/fix_nve_kokkos.h +++ b/src/KOKKOS/fix_nve_kokkos.h @@ -46,6 +46,7 @@ class FixNVEKokkos : public FixNVE { void init() override; void initial_integrate(int) override; void final_integrate() override; + void fused_integrate() override; KOKKOS_INLINE_FUNCTION void initial_integrate_item(int) const; @@ -96,6 +97,25 @@ struct FixNVEKokkosFinalIntegrateFunctor { } }; +template +struct FixNVEKokkosFusedIntegrateFunctor { + typedef DeviceType device_type ; + FixNVEKokkos c; + + FixNVEKokkosFusedIntegrateFunctor(FixNVEKokkos* c_ptr): + c(*c_ptr) {c.cleanup_copy();}; + KOKKOS_INLINE_FUNCTION + void operator()(const int i) const { + if (RMass) { + c.final_integrate_rmass_item(i); + c.initial_integrate_rmass_item(i); + } else { + c.final_integrate_item(i); + c.initial_integrate_item(i); + } + } +}; + } #endif diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp index 9d8c16603e..581f6598f1 100644 --- a/src/KOKKOS/modify_kokkos.cpp +++ b/src/KOKKOS/modify_kokkos.cpp @@ -392,6 +392,25 @@ void ModifyKokkos::final_integrate() } } + +/* ---------------------------------------------------------------------- + 2nd half of integrate call, only for relevant fixes +------------------------------------------------------------------------- */ + +void ModifyKokkos::fused_integrate() +{ + for (int i = 0; i < n_final_integrate; i++) { + atomKK->sync(fix[list_final_integrate[i]]->execution_space, + fix[list_final_integrate[i]]->datamask_read); + int prev_auto_sync = lmp->kokkos->auto_sync; + if (!fix[list_final_integrate[i]]->kokkosable) lmp->kokkos->auto_sync = 1; + fix[list_final_integrate[i]]->fused_integrate(); + lmp->kokkos->auto_sync = prev_auto_sync; + atomKK->modified(fix[list_final_integrate[i]]->execution_space, + fix[list_final_integrate[i]]->datamask_modify); + } +} + /* ---------------------------------------------------------------------- end-of-timestep call, only for relevant fixes only call fix->end_of_step() on timesteps that are multiples of nevery diff --git a/src/KOKKOS/modify_kokkos.h b/src/KOKKOS/modify_kokkos.h index 5edf5cd662..e440693fbb 100644 --- a/src/KOKKOS/modify_kokkos.h +++ b/src/KOKKOS/modify_kokkos.h @@ -39,6 +39,7 @@ class ModifyKokkos : public Modify { void pre_reverse(int,int) override; void post_force(int) override; void final_integrate() override; + void fused_integrate() override; void end_of_step() override; double energy_couple() override; double energy_global() override; diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h index 0ff244f67d..7551d03f89 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/src/KOKKOS/pair_kokkos.h @@ -280,6 +280,12 @@ struct PairComputeFunctor { const X_FLOAT ztmp = c.x(i,2); const int itype = c.type(i); + Kokkos::single(Kokkos::PerThread(team), [&] (){ + f(i,0) = 0.0; + f(i,1) = 0.0; + f(i,2) = 0.0; + }); + const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i); const int jnum = list.d_numneigh[i]; @@ -337,6 +343,12 @@ struct PairComputeFunctor { const int itype = c.type(i); const F_FLOAT qtmp = c.q(i); + Kokkos::single(Kokkos::PerThread(team), [&] (){ + f(i,0) = 0.0; + f(i,1) = 0.0; + f(i,2) = 0.0; + }); + const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i); const int jnum = list.d_numneigh[i]; @@ -399,6 +411,12 @@ struct PairComputeFunctor { const X_FLOAT ztmp = c.x(i,2); const int itype = c.type(i); + Kokkos::single(Kokkos::PerThread(team), [&] (){ + f(i,0) = 0.0; + f(i,1) = 0.0; + f(i,2) = 0.0; + }); + const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i); const int jnum = list.d_numneigh[i]; @@ -495,6 +513,12 @@ struct PairComputeFunctor { const int itype = c.type(i); const F_FLOAT qtmp = c.q(i); + Kokkos::single(Kokkos::PerThread(team), [&] (){ + f(i,0) = 0.0; + f(i,1) = 0.0; + f(i,2) = 0.0; + }); + const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i); const int jnum = list.d_numneigh[i]; @@ -743,6 +767,8 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename std::enable_if<(NEIG fpair->lmp->kokkos->neigh_thread = 1; if (fpair->lmp->kokkos->neigh_thread) { + fpair->fuse_force_clear_flag = 1; + int vector_length = 8; int atoms_per_team = 32; diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index 01401542cd..0339df663b 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -296,7 +296,9 @@ void VerletKokkos::run(int n) // initial time integration timer->stamp(); - modify->initial_integrate(vflag); + fuse_check(i,n); + if (!fuse_integrate) + modify->initial_integrate(vflag); if (n_post_integrate) modify->post_integrate(); timer->stamp(Timer::MODIFY); @@ -362,7 +364,8 @@ void VerletKokkos::run(int n) // since some bonded potentials tally pairwise energy/virial // and Pair:ev_tally() needs to be called before any tallying - force_clear(); + if (!fuse_force_clear) + force_clear(); timer->stamp(); @@ -494,7 +497,10 @@ void VerletKokkos::run(int n) // force modifications, final time integration, diagnostics if (n_post_force) modify->post_force(vflag); - modify->final_integrate(); + + if (fuse_integrate) modify->fused_integrate(); + else modify->final_integrate(); + if (n_end_of_step) modify->end_of_step(); timer->stamp(Timer::MODIFY); @@ -593,3 +599,35 @@ void VerletKokkos::force_clear() } } } + +/* ---------------------------------------------------------------------- + check if can fuse force_clear() with pair compute() + Requirements: + - no pre_force fixes + - no torques, SPIN forces, or includegroup set + - pair compute() must be called + - pair_style must support fusing + + check if can fuse initial_integrate() with final_integrate() + Requirements: + - no end_of_step fixes + - not on first, last, or output step + - no timers to break out of loop + - integrate fix style must support fusing +------------------------------------------------------------------------- */ + +void VerletKokkos::fuse_check(int i, int n) +{ + fuse_force_clear = 0; + if (modify->n_pre_force) fuse_force_clear = 0; + if (torqueflag || extraflag || neighbor->includegroup) fuse_force_clear = 0; + if (!pair_compute_flag) fuse_force_clear = 0; + if (!force->pair->fuse_force_clear_flag) fuse_force_clear = 0; + + fuse_integrate = 0; + if (modify->n_end_of_step) fuse_integrate = 0; + if (i == 0 || i == n-1) fuse_integrate = 0; + if (update->ntimestep == output->next) fuse_integrate = 0; + if (timer->has_timeout()) fuse_integrate = 0; + if (!modify->check_fuse_integrate()) fuse_integrate = 0; +} diff --git a/src/KOKKOS/verlet_kokkos.h b/src/KOKKOS/verlet_kokkos.h index 067df54f4f..c71211c542 100644 --- a/src/KOKKOS/verlet_kokkos.h +++ b/src/KOKKOS/verlet_kokkos.h @@ -46,6 +46,9 @@ class VerletKokkos : public Verlet { protected: DAT::t_f_array f_merge_copy,f; + int fuse_force_clear,fuse_integrate; + + void fuse_check(int, int); }; } diff --git a/src/fix.cpp b/src/fix.cpp index 1d41ad3943..02adcbd016 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -102,15 +102,15 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) : vflag_atom = cvflag_atom = 0; centroidstressflag = CENTROID_SAME; - // KOKKOS per-fix data masks + // KOKKOS package execution_space = Host; datamask_read = ALL_MASK; datamask_modify = ALL_MASK; - kokkosable = 0; + kokkosable = copymode = 0; forward_comm_device = exchange_comm_device = 0; - copymode = 0; + fuse_integrate_flag = 0; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix.h b/src/fix.h index b47cfb2f4a..23b120d989 100644 --- a/src/fix.h +++ b/src/fix.h @@ -127,11 +127,12 @@ class Fix : protected Pointers { int restart_reset; // 1 if restart just re-initialized fix - // KOKKOS host/device flag and data masks + // KOKKOS flags and variables int kokkosable; // 1 if Kokkos fix int forward_comm_device; // 1 if forward comm on Device int exchange_comm_device; // 1 if exchange comm on Device + int fuse_integrate_flag; // 1 if can fuse initial integrate with final integrate ExecutionSpace execution_space; unsigned int datamask_read, datamask_modify; @@ -152,6 +153,7 @@ class Fix : protected Pointers { virtual void setup_pre_reverse(int, int) {} virtual void min_setup(int) {} virtual void initial_integrate(int) {} + virtual void fused_integrate() {} virtual void post_integrate() {} virtual void pre_exchange() {} virtual void pre_neighbor() {} diff --git a/src/modify.cpp b/src/modify.cpp index d0656d3895..71f7cb8889 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -475,6 +475,17 @@ void Modify::final_integrate() for (int i = 0; i < n_final_integrate; i++) fix[list_final_integrate[i]]->final_integrate(); } + +/* ---------------------------------------------------------------------- + 2nd half of integrate call, only for relevant fixes +------------------------------------------------------------------------- */ + +void Modify::fused_integrate() +{ + for (int i = 0; i < n_final_integrate; i++) + fix[list_final_integrate[i]]->fused_integrate(); +} + /* ---------------------------------------------------------------------- end-of-timestep call, only for relevant fixes only call fix->end_of_step() on timesteps that are multiples of nevery @@ -1799,3 +1810,22 @@ double Modify::memory_usage() for (int i = 0; i < ncompute; i++) bytes += compute[i]->memory_usage(); return bytes; } + +/* ---------------------------------------------------------------------- + check if initial and final integrate can be fused +------------------------------------------------------------------------- */ + +int Modify::check_fuse_integrate() +{ + int fuse_integrate_flag = 1; + + for (int i = 0; i < n_initial_integrate; i++) + if (!fix[list_initial_integrate[i]]->fuse_integrate_flag) + fuse_integrate_flag = 0; + + for (int i = 0; i < n_final_integrate; i++) + if (!fix[list_final_integrate[i]]->fuse_integrate_flag) + fuse_integrate_flag = 0; + + return fuse_integrate_flag; +} diff --git a/src/modify.h b/src/modify.h index 7a3f54c277..da3cab55d2 100644 --- a/src/modify.h +++ b/src/modify.h @@ -61,6 +61,7 @@ class Modify : protected Pointers { virtual void setup_pre_force(int); virtual void setup_pre_reverse(int, int); virtual void initial_integrate(int); + virtual void fused_integrate(); virtual void post_integrate(); virtual void pre_exchange(); virtual void pre_neighbor(); @@ -150,6 +151,8 @@ class Modify : protected Pointers { double memory_usage(); + int check_fuse_integrate(); + protected: // internal fix counts diff --git a/src/pair.cpp b/src/pair.cpp index 9ae24e1e93..34c8e4978e 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -116,15 +116,14 @@ Pair::Pair(LAMMPS *lmp) : nondefault_history_transfer = 0; beyond_contact = 0; - // KOKKOS per-fix data masks + // KOKKOS package execution_space = Host; datamask_read = ALL_MASK; datamask_modify = ALL_MASK; - kokkosable = 0; - reverse_comm_device = 0; - copymode = 0; + kokkosable = copymode = 0; + reverse_comm_device = fuse_force_clear_flag = 0; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair.h b/src/pair.h index 8c856660e9..d36b48f340 100644 --- a/src/pair.h +++ b/src/pair.h @@ -119,12 +119,13 @@ class Pair : protected Pointers { int beyond_contact, nondefault_history_transfer; // for granular styles - // KOKKOS host/device flag and data masks + // KOKKOS flags and variables ExecutionSpace execution_space; unsigned int datamask_read, datamask_modify; int kokkosable; // 1 if Kokkos pair int reverse_comm_device; // 1 if reverse comm on Device + int fuse_force_clear_flag; // 1 if can fuse force clear with force compute Pair(class LAMMPS *); ~Pair() override; diff --git a/src/timer.h b/src/timer.h index 5c100db1c0..f7efa5ac64 100644 --- a/src/timer.h +++ b/src/timer.h @@ -63,6 +63,7 @@ class Timer : protected Pointers { bool has_normal() const { return (_level >= NORMAL); } bool has_full() const { return (_level >= FULL); } bool has_sync() const { return (_sync != OFF); } + bool has_timeout() const { return (_timeout >= 0.0); } // flag if wallclock time is expired bool is_timeout() const { return (_timeout == 0.0); } From 29a2b58ab625b9ea6361a0e82a14369ec91fbc68 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Fri, 28 Apr 2023 14:57:19 -0600 Subject: [PATCH 02/41] Small cleanup --- src/KOKKOS/modify_kokkos.cpp | 21 ++++++++++++++++++++- src/KOKKOS/modify_kokkos.h | 2 ++ src/KOKKOS/verlet_kokkos.cpp | 4 ++-- src/modify.cpp | 30 ------------------------------ src/modify.h | 4 +--- 5 files changed, 25 insertions(+), 36 deletions(-) diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp index 581f6598f1..d0f57df8dc 100644 --- a/src/KOKKOS/modify_kokkos.cpp +++ b/src/KOKKOS/modify_kokkos.cpp @@ -394,7 +394,7 @@ void ModifyKokkos::final_integrate() /* ---------------------------------------------------------------------- - 2nd half of integrate call, only for relevant fixes + fused initial and final integrate call, only for relevant fixes ------------------------------------------------------------------------- */ void ModifyKokkos::fused_integrate() @@ -900,3 +900,22 @@ int ModifyKokkos::min_reset_ref() } return itmpall; } + +/* ---------------------------------------------------------------------- + check if initial and final integrate can be fused +------------------------------------------------------------------------- */ + +int ModifyKokkos::check_fuse_integrate() +{ + int fuse_integrate_flag = 1; + + for (int i = 0; i < n_initial_integrate; i++) + if (!fix[list_initial_integrate[i]]->fuse_integrate_flag) + fuse_integrate_flag = 0; + + for (int i = 0; i < n_final_integrate; i++) + if (!fix[list_final_integrate[i]]->fuse_integrate_flag) + fuse_integrate_flag = 0; + + return fuse_integrate_flag; +} diff --git a/src/KOKKOS/modify_kokkos.h b/src/KOKKOS/modify_kokkos.h index e440693fbb..8dd9e6d9df 100644 --- a/src/KOKKOS/modify_kokkos.h +++ b/src/KOKKOS/modify_kokkos.h @@ -70,6 +70,8 @@ class ModifyKokkos : public Modify { int min_dof() override; int min_reset_ref() override; + int check_fuse_integrate(); + protected: }; diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index 0339df663b..2b70f3db61 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -27,7 +27,7 @@ #include "kspace.h" #include "output.h" #include "update.h" -#include "modify.h" +#include "modify_kokkos.h" #include "timer.h" #include "memory_kokkos.h" #include "kokkos.h" @@ -629,5 +629,5 @@ void VerletKokkos::fuse_check(int i, int n) if (i == 0 || i == n-1) fuse_integrate = 0; if (update->ntimestep == output->next) fuse_integrate = 0; if (timer->has_timeout()) fuse_integrate = 0; - if (!modify->check_fuse_integrate()) fuse_integrate = 0; + if (!((ModifyKokkos*)modify)->check_fuse_integrate()) fuse_integrate = 0; } diff --git a/src/modify.cpp b/src/modify.cpp index 71f7cb8889..d0656d3895 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -475,17 +475,6 @@ void Modify::final_integrate() for (int i = 0; i < n_final_integrate; i++) fix[list_final_integrate[i]]->final_integrate(); } - -/* ---------------------------------------------------------------------- - 2nd half of integrate call, only for relevant fixes -------------------------------------------------------------------------- */ - -void Modify::fused_integrate() -{ - for (int i = 0; i < n_final_integrate; i++) - fix[list_final_integrate[i]]->fused_integrate(); -} - /* ---------------------------------------------------------------------- end-of-timestep call, only for relevant fixes only call fix->end_of_step() on timesteps that are multiples of nevery @@ -1810,22 +1799,3 @@ double Modify::memory_usage() for (int i = 0; i < ncompute; i++) bytes += compute[i]->memory_usage(); return bytes; } - -/* ---------------------------------------------------------------------- - check if initial and final integrate can be fused -------------------------------------------------------------------------- */ - -int Modify::check_fuse_integrate() -{ - int fuse_integrate_flag = 1; - - for (int i = 0; i < n_initial_integrate; i++) - if (!fix[list_initial_integrate[i]]->fuse_integrate_flag) - fuse_integrate_flag = 0; - - for (int i = 0; i < n_final_integrate; i++) - if (!fix[list_final_integrate[i]]->fuse_integrate_flag) - fuse_integrate_flag = 0; - - return fuse_integrate_flag; -} diff --git a/src/modify.h b/src/modify.h index da3cab55d2..9686a9a7ec 100644 --- a/src/modify.h +++ b/src/modify.h @@ -61,7 +61,7 @@ class Modify : protected Pointers { virtual void setup_pre_force(int); virtual void setup_pre_reverse(int, int); virtual void initial_integrate(int); - virtual void fused_integrate(); + virtual void fused_integrate() {} virtual void post_integrate(); virtual void pre_exchange(); virtual void pre_neighbor(); @@ -151,8 +151,6 @@ class Modify : protected Pointers { double memory_usage(); - int check_fuse_integrate(); - protected: // internal fix counts From 24e7204b51c7392e1dc556394028b76d163354a5 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Fri, 28 Apr 2023 15:00:50 -0600 Subject: [PATCH 03/41] whitespace --- src/KOKKOS/modify_kokkos.cpp | 1 - src/KOKKOS/pair_kokkos.h | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp index d0f57df8dc..44fcb48727 100644 --- a/src/KOKKOS/modify_kokkos.cpp +++ b/src/KOKKOS/modify_kokkos.cpp @@ -392,7 +392,6 @@ void ModifyKokkos::final_integrate() } } - /* ---------------------------------------------------------------------- fused initial and final integrate call, only for relevant fixes ------------------------------------------------------------------------- */ diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h index 7551d03f89..dd274ca8c2 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/src/KOKKOS/pair_kokkos.h @@ -416,7 +416,7 @@ struct PairComputeFunctor { f(i,1) = 0.0; f(i,2) = 0.0; }); - + const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i); const int jnum = list.d_numneigh[i]; From f66cf7c9e8a8a02d552e7a7dd489d1acc1f3a334 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Fri, 28 Apr 2023 16:15:55 -0600 Subject: [PATCH 04/41] Revert testing change --- src/KOKKOS/verlet_kokkos.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index 2b70f3db61..4e564a6ef2 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -618,13 +618,13 @@ void VerletKokkos::force_clear() void VerletKokkos::fuse_check(int i, int n) { - fuse_force_clear = 0; + fuse_force_clear = 1; if (modify->n_pre_force) fuse_force_clear = 0; if (torqueflag || extraflag || neighbor->includegroup) fuse_force_clear = 0; if (!pair_compute_flag) fuse_force_clear = 0; if (!force->pair->fuse_force_clear_flag) fuse_force_clear = 0; - fuse_integrate = 0; + fuse_integrate = 1; if (modify->n_end_of_step) fuse_integrate = 0; if (i == 0 || i == n-1) fuse_integrate = 0; if (update->ntimestep == output->next) fuse_integrate = 0; From 2c8a49bb2600d4ebcf54c8a991a72d2845d36d8a Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Sat, 29 Apr 2023 18:14:24 -0600 Subject: [PATCH 05/41] new compute count_type --- src/compute_count_type.cpp | 165 +++++++++++++++++++++++++++++++++++++ src/compute_count_type.h | 46 +++++++++++ src/lmptype.h | 5 +- 3 files changed, 215 insertions(+), 1 deletion(-) create mode 100644 src/compute_count_type.cpp create mode 100644 src/compute_count_type.h diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp new file mode 100644 index 0000000000..707a4bb943 --- /dev/null +++ b/src/compute_count_type.cpp @@ -0,0 +1,165 @@ +/* ---------------------------------------------------------------------- + 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 "compute_count_type.h" + +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "group.h" +#include "update.h" + +using namespace LAMMPS_NS; + +enum{ATOM,BOND}; + +/* ---------------------------------------------------------------------- */ + +ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR, "Illegal compute count/type command"); + + // process args + + if (strcmp(arg[3],"atom") == 0) mode = ATOM; + else if (strcmp(arg[3],"bond") == 0) mode = BOND; + else error->all(FLERR, "Invalid compute count/type keyword {}",arg[3]); + + if (mode == ATOM) { + vector_flag = 1; + size_vector = atom->ntypes; + extvector = 1; + } else if (mode == BOND) { + scalar_flag = vector_flag = 1; + size_vector = atom->nbondtypes; + extscalar = 1; + extvector = 1; + } + + if (mode == BOND && !atom->avec->bonds_allow) + error->all(FLERR,"Cannot use compute count/type bond command with no bonds allowed"); + + vector = new double[size_vector]; + + // work vectors + + count = new int[size_vector]; + bcount_me = new bigint[size_vector]; + bcount = new bigint[size_vector]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeCountType::~ComputeCountType() +{ + delete[] vector; +} + +/* ---------------------------------------------------------------------- + only invoked for mode = BOND to count broken bonds + broken bonds have bond_type = 0 +---------------------------------------------------------------------- */ + +double ComputeCountType::compute_scalar() +{ + invoked_scalar = update->ntimestep; + + int *num_bond = atom->num_bond; + int **bond_type = atom->bond_type; + int nlocal = atom->nlocal; + + int nbond; + int count = 0; + + for (int i = 0; i < nlocal; i++) { + nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) + if (bond_type[i][m] == 0) count++; + } + + // sum across procs as bigint, then convert to double + // correct for double counting if newton_bond off + + bigint bcount_me = count; + bigint bcount; + MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); + if (force->newton_bond == 0) bcount /= 2; + + if (bcount > MAXDOUBLEINT) + error->all(FLERR,"Compute count/type overflow"); + double scalar = bcount; + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCountType::compute_vector() +{ + invoked_vector = update->ntimestep; + + int n; + + // count atoms by type + + if (mode == ATOM) { + int *type = atom->type; + int nlocal = atom->nlocal; + int ntypes = atom->ntypes; + + for (int m = 0; m < ntypes; m++) count[m] = 0; + for (int i = 0; i < nlocal; i++) count[type[i]-1]++; + + n = ntypes; + } + + // count bonds by type + // skip type = 0 bonds, they are counted by compute_scalar + // bond types can be negative for SHAKE + + else if (mode == BOND) { + int *num_bond = atom->num_bond; + int **bond_type = atom->bond_type; + int nlocal = atom->nlocal; + int nbondtypes = atom->nbondtypes; + + int nbond,itype; + for (int m = 0; m < nbondtypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) { + itype = bond_type[i][m]; + if (itype == 0) continue; + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + + n = nbondtypes; + } + + // sum across procs as bigint, then convert to double + // correct for double counting if newton_bond off + + for (int m = 0; m < n; m++) bcount_me[m] = count[m]; + MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); + if (force->newton_bond == 0) + for (int m = 0; m < n; m++) bcount[m] /= 2; + + for (int m = 0; m < n; m++) + if (bcount[m] > MAXDOUBLEINT) + error->all(FLERR,"Compute count/type overflow"); + for (int m = 0; m < n; m++) vector[m] *= bcount[m]; +} diff --git a/src/compute_count_type.h b/src/compute_count_type.h new file mode 100644 index 0000000000..0ed7a22467 --- /dev/null +++ b/src/compute_count_type.h @@ -0,0 +1,46 @@ +/* -*- 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 COMPUTE_CLASS +// clang-format off +ComputeStyle(count/type,ComputeCountType); +// clang-format on +#else + +#ifndef LMP_COMPUTE_COMPUTE_TYPE_H +#define LMP_COMPUTE_COMPUTE_TYPE_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeCountType : public Compute { + public: + ComputeCountType(class LAMMPS *, int, char **); + ~ComputeCountType() override; + void init() override {} + double compute_scalar() override; + void compute_vector() override; + + protected: + int mode; + + int *count; + bigint *bcount_me; + bigint *bcount; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/lmptype.h b/src/lmptype.h index bf56d07d77..891d5bdf89 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -95,7 +95,8 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT_MAX #define MAXBIGINT INT64_MAX - +#define MAXDOUBLEINT 9007199254740992 // 2^53 + #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT #define MPI_LMP_BIGINT MPI_LL @@ -132,6 +133,7 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT64_MAX #define MAXBIGINT INT64_MAX +#define MAXDOUBLEINT 9007199254740992 // 2^53 #define MPI_LMP_TAGINT MPI_LL #define MPI_LMP_IMAGEINT MPI_LL @@ -168,6 +170,7 @@ typedef int bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT_MAX #define MAXBIGINT INT_MAX +#define MAXDOUBLEINT INT_MAX #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT From f644f7078cb782702df393dafce4248c36c84ff9 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Sat, 29 Apr 2023 18:50:15 -0600 Subject: [PATCH 06/41] doc page for new command --- doc/src/Commands_compute.rst | 1 + doc/src/bond_style.rst | 14 +++--- doc/src/compute.rst | 1 + doc/src/compute_count_type.rst | 89 ++++++++++++++++++++++++++++++++++ src/compute_count_type.cpp | 8 ++- 5 files changed, 105 insertions(+), 8 deletions(-) create mode 100644 doc/src/compute_count_type.rst diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 2f9c8b03fb..755000c976 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -46,6 +46,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`com/chunk ` * :doc:`contact/atom ` * :doc:`coord/atom (k) ` + * :doc:`count/type ` * :doc:`damage/atom ` * :doc:`dihedral ` * :doc:`dihedral/local ` diff --git a/doc/src/bond_style.rst b/doc/src/bond_style.rst index b33d0a9e9a..95f463e695 100644 --- a/doc/src/bond_style.rst +++ b/doc/src/bond_style.rst @@ -32,13 +32,13 @@ Set the formula(s) LAMMPS uses to compute bond interactions between pairs of atoms. In LAMMPS, a bond differs from a pairwise interaction, which are set via the :doc:`pair_style ` command. Bonds are defined between specified pairs of atoms and -remain in force for the duration of the simulation (unless the bond -breaks which is possible in some bond potentials). The list of bonded -atoms is read in by a :doc:`read_data ` or -:doc:`read_restart ` command from a data or restart file. -By contrast, pair potentials are typically defined between all pairs -of atoms within a cutoff distance and the set of active interactions -changes over time. +remain in force for the duration of the simulation (unless new bonds +are created or existing bonds break, which is possible in some fixes +and bond potentials). The list of bonded atoms is read in by a +:doc:`read_data ` or :doc:`read_restart ` +command from a data or restart file. By contrast, pair potentials are +typically defined between all pairs of atoms within a cutoff distance +and the set of active interactions changes over time. Hybrid models where bonds are computed using different bond potentials can be setup using the *hybrid* bond style. diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 880f60a8a6..950487cb72 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -200,6 +200,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`com/chunk ` - center of mass for each chunk * :doc:`contact/atom ` - contact count for each spherical particle * :doc:`coord/atom ` - coordination number for each atom +* :doc:`count/type ` - count of atoms or bonds by type * :doc:`damage/atom ` - Peridynamic damage for each atom * :doc:`dihedral ` - energy of each dihedral sub-style * :doc:`dihedral/local ` - angle of each dihedral diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst new file mode 100644 index 0000000000..17fd9ed784 --- /dev/null +++ b/doc/src/compute_count_type.rst @@ -0,0 +1,89 @@ +.. index:: compute count/type + +compute count/type command +==================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + compute ID group-ID count/type mode + +* ID, group-ID are documented in :doc:`compute ` command +* count/type = style name of this compute command +* mode = {atom} or {bond} + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all count/type atom + compute 1 flowmols count/type bond + +Description +""""""""""" + +Define a computation that counts the current number of atoms by atom +type or the number of bonds by bond type. The latter can be useful in +reactive simulations where bonds are broken or created. + +Note that for this command, bonds are the topological kind enumerated +in a data file, initially read by the :doc:`read_data ` +command. They do not refer to bonds defined on-the-fly by bond-order +or reactive pair styles. + +These commands can create and break toplogical bonds: + +* :doc:`fix bond/react ` +* :doc:`fix bond/create ` +* :doc:`fix bond/break ` +* :doc:`bond_style quartic ` +* :doc:`BPM package ` bond styles + +If the {mode} setting is {atom} then the count of atoms for each atom +type is tallied. Only atoms in the specified group are counted. + +If the {mode} setting is {bond} then the count of bonds for each bond +type is tallied. Only bonds with both atoms in the specified group +are counted. + +For {mode} = {bond}, broken bonds with a bond type or zero are also +counted. Some commands flag broken bonds by setting their bond type +to zero. See the :doc:`Howto broken bonds ` doc +page for details. Note that the group setting is ignored for broken +bonds; all broken bonds in the system are counted. + +---------- + +Output info +""""""""""" + +This compute calculates a global vector of counts. If the mode is +{atom}, the vector length is the number of atom types. If the mode is +{bond}, the vector length is the number of bond types. + +If the mode is {bond} this compute also calculates a global scalar which +counts the number of broken bonds. + +These values can be used by any command that uses global scalar or +vector values from a compute as input. See the :doc:`Howto output +` page for an overview of LAMMPS output options. + +The scalar and vector values calculated by this compute are "extensive". + +Restrictions +"""""""""""" + +none + +Related commands +"""""""""""""""" + +none + +Default +""""""" + +none diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 707a4bb943..c3014394d8 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -82,6 +82,8 @@ double ComputeCountType::compute_scalar() int nbond; int count = 0; + + // NOTE: respect group setting for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; @@ -113,6 +115,8 @@ void ComputeCountType::compute_vector() // count atoms by type + // NOTE: respect group setting + if (mode == ATOM) { int *type = atom->type; int nlocal = atom->nlocal; @@ -127,7 +131,9 @@ void ComputeCountType::compute_vector() // count bonds by type // skip type = 0 bonds, they are counted by compute_scalar // bond types can be negative for SHAKE - + + // NOTE: respect group setting + else if (mode == BOND) { int *num_bond = atom->num_bond; int **bond_type = atom->bond_type; From c282d8d5d0eb80c9050f8b141639ef1e31617e06 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 1 May 2023 09:01:03 -0600 Subject: [PATCH 07/41] add group support to new command --- doc/src/compute_count_type.rst | 11 ++++----- src/compute_count_type.cpp | 45 ++++++++++++++++++++++++---------- 2 files changed, 37 insertions(+), 19 deletions(-) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index 17fd9ed784..ca56648297 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -29,7 +29,7 @@ Define a computation that counts the current number of atoms by atom type or the number of bonds by bond type. The latter can be useful in reactive simulations where bonds are broken or created. -Note that for this command, bonds are the topological kind enumerated +Note that for this command, bonds are the topological ones enumerated in a data file, initially read by the :doc:`read_data ` command. They do not refer to bonds defined on-the-fly by bond-order or reactive pair styles. @@ -49,9 +49,8 @@ If the {mode} setting is {bond} then the count of bonds for each bond type is tallied. Only bonds with both atoms in the specified group are counted. -For {mode} = {bond}, broken bonds with a bond type or zero are also -counted. Some commands flag broken bonds by setting their bond type -to zero. See the :doc:`Howto broken bonds ` doc +For {mode} = {bond}, broken bonds with a bond type of zero are also +counted. See the :doc:`Howto broken bonds ` doc page for details. Note that the group setting is ignored for broken bonds; all broken bonds in the system are counted. @@ -64,8 +63,8 @@ This compute calculates a global vector of counts. If the mode is {atom}, the vector length is the number of atom types. If the mode is {bond}, the vector length is the number of bond types. -If the mode is {bond} this compute also calculates a global scalar which -counts the number of broken bonds. +If the mode is {bond} this compute also calculates a global scalar +which is the number of broken bonds. These values can be used by any command that uses global scalar or vector values from a compute as input. See the :doc:`Howto output diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index c3014394d8..9de5f39cab 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -80,14 +80,16 @@ double ComputeCountType::compute_scalar() int **bond_type = atom->bond_type; int nlocal = atom->nlocal; - int nbond; + int m,nbond; int count = 0; - // NOTE: respect group setting + // count broken bonds with bond_type = 0 + // ignore group setting since 2 atoms in a broken bond + // can be arbitrarily far apart for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; - for (int m = 0; m < nbond; m++) + for (m = 0; m < nbond; m++) if (bond_type[i][m] == 0) count++; } @@ -114,45 +116,62 @@ void ComputeCountType::compute_vector() int n; // count atoms by type - - // NOTE: respect group setting + // atom must be in group to be counted if (mode == ATOM) { int *type = atom->type; + int *mask = atom->mask; int nlocal = atom->nlocal; int ntypes = atom->ntypes; for (int m = 0; m < ntypes; m++) count[m] = 0; - for (int i = 0; i < nlocal; i++) count[type[i]-1]++; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + count[type[i]-1]++; n = ntypes; } // count bonds by type + // both atoms in bond must be in group to be counted // skip type = 0 bonds, they are counted by compute_scalar // bond types can be negative for SHAKE - // NOTE: respect group setting - else if (mode == BOND) { - int *num_bond = atom->num_bond; + int **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; + int *num_bond = atom->num_bond; + int *mask = atom->mask; int nlocal = atom->nlocal; int nbondtypes = atom->nbondtypes; - int nbond,itype; + int j,m,nbond,itype; + int flag = 0; for (int m = 0; m < nbondtypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; - for (int m = 0; m < nbond; m++) { + for (m = 0; m < nbond; m++) { itype = bond_type[i][m]; if (itype == 0) continue; - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + + j = atom->map(bond_atom[i][m]); + if (j < 0) { + flag = 1; + continue; + } + + if ((mask[i] & groupbit) && (mask[j] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } } } + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); + n = nbondtypes; } From 2fe423cc70a104514b8f2c396e3b732a86c74d0d Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 1 May 2023 09:20:20 -0600 Subject: [PATCH 08/41] bug fix --- src/compute_count_type.cpp | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 9de5f39cab..feadcb04ba 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -51,6 +51,8 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( if (mode == BOND && !atom->avec->bonds_allow) error->all(FLERR,"Cannot use compute count/type bond command with no bonds allowed"); + // output vector + vector = new double[size_vector]; // work vectors @@ -65,6 +67,10 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( ComputeCountType::~ComputeCountType() { delete[] vector; + + delete[] count; + delete[] bcount_me; + delete[] bcount; } /* ---------------------------------------------------------------------- @@ -103,7 +109,7 @@ double ComputeCountType::compute_scalar() if (bcount > MAXDOUBLEINT) error->all(FLERR,"Compute count/type overflow"); - double scalar = bcount; + scalar = bcount; return scalar; } @@ -113,7 +119,7 @@ void ComputeCountType::compute_vector() { invoked_vector = update->ntimestep; - int n; + int nvec; // count atoms by type // atom must be in group to be counted @@ -129,7 +135,7 @@ void ComputeCountType::compute_vector() if (mask[i] & groupbit) count[type[i]-1]++; - n = ntypes; + nvec = ntypes; } // count bonds by type @@ -172,19 +178,19 @@ void ComputeCountType::compute_vector() MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); - n = nbondtypes; + nvec = nbondtypes; } // sum across procs as bigint, then convert to double // correct for double counting if newton_bond off - for (int m = 0; m < n; m++) bcount_me[m] = count[m]; - MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); + for (int m = 0; m < nvec; m++) bcount_me[m] = count[m]; + MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world); if (force->newton_bond == 0) - for (int m = 0; m < n; m++) bcount[m] /= 2; + for (int m = 0; m < nvec; m++) bcount[m] /= 2; - for (int m = 0; m < n; m++) + for (int m = 0; m < nvec; m++) if (bcount[m] > MAXDOUBLEINT) error->all(FLERR,"Compute count/type overflow"); - for (int m = 0; m < n; m++) vector[m] *= bcount[m]; + for (int m = 0; m < nvec; m++) vector[m] = bcount[m]; } From 370ec562072ca886df77f8978dea407809d5fa7f Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 1 May 2023 10:20:45 -0600 Subject: [PATCH 09/41] Whitespace --- src/KOKKOS/pair_kokkos.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h index dd274ca8c2..c244d6a8f9 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/src/KOKKOS/pair_kokkos.h @@ -416,7 +416,7 @@ struct PairComputeFunctor { f(i,1) = 0.0; f(i,2) = 0.0; }); - + const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i); const int jnum = list.d_numneigh[i]; From 059601734f8fee570af87fff1d961212be1d3299 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 1 May 2023 12:38:16 -0600 Subject: [PATCH 10/41] Fix issues --- src/KOKKOS/fix_nve_kokkos.cpp | 3 ++- src/KOKKOS/verlet_kokkos.cpp | 18 +++++++++--------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/KOKKOS/fix_nve_kokkos.cpp b/src/KOKKOS/fix_nve_kokkos.cpp index f812c575bc..c26e26a02c 100644 --- a/src/KOKKOS/fix_nve_kokkos.cpp +++ b/src/KOKKOS/fix_nve_kokkos.cpp @@ -177,7 +177,6 @@ template void FixNVEKokkos::fused_integrate() { atomKK->sync(execution_space,datamask_read); - atomKK->modified(execution_space,datamask_modify); x = atomKK->k_x.view(); v = atomKK->k_v.view(); @@ -196,6 +195,8 @@ void FixNVEKokkos::fused_integrate() FixNVEKokkosFusedIntegrateFunctor functor(this); Kokkos::parallel_for(nlocal,functor); } + + atomKK->modified(execution_space,datamask_modify); } namespace LAMMPS_NS { diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index 4e564a6ef2..bd15859644 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -296,7 +296,6 @@ void VerletKokkos::run(int n) // initial time integration timer->stamp(); - fuse_check(i,n); if (!fuse_integrate) modify->initial_integrate(vflag); if (n_post_integrate) modify->post_integrate(); @@ -364,6 +363,7 @@ void VerletKokkos::run(int n) // since some bonded potentials tally pairwise energy/virial // and Pair:ev_tally() needs to be called before any tallying + fuse_check(i,n); if (!fuse_force_clear) force_clear(); @@ -611,7 +611,7 @@ void VerletKokkos::force_clear() check if can fuse initial_integrate() with final_integrate() Requirements: - no end_of_step fixes - - not on first, last, or output step + - not on last or output step - no timers to break out of loop - integrate fix style must support fusing ------------------------------------------------------------------------- */ @@ -620,14 +620,14 @@ void VerletKokkos::fuse_check(int i, int n) { fuse_force_clear = 1; if (modify->n_pre_force) fuse_force_clear = 0; - if (torqueflag || extraflag || neighbor->includegroup) fuse_force_clear = 0; - if (!pair_compute_flag) fuse_force_clear = 0; - if (!force->pair->fuse_force_clear_flag) fuse_force_clear = 0; + else if (torqueflag || extraflag || neighbor->includegroup) fuse_force_clear = 0; + else if (!force->pair || !pair_compute_flag) fuse_force_clear = 0; + else if (!force->pair->fuse_force_clear_flag) fuse_force_clear = 0; fuse_integrate = 1; if (modify->n_end_of_step) fuse_integrate = 0; - if (i == 0 || i == n-1) fuse_integrate = 0; - if (update->ntimestep == output->next) fuse_integrate = 0; - if (timer->has_timeout()) fuse_integrate = 0; - if (!((ModifyKokkos*)modify)->check_fuse_integrate()) fuse_integrate = 0; + else if (i == n-1) fuse_integrate = 0; + else if (update->ntimestep == output->next) fuse_integrate = 0; + else if (timer->has_timeout()) fuse_integrate = 0; + else if (!((ModifyKokkos*)modify)->check_fuse_integrate()) fuse_integrate = 0; } From fd65000afc950d7cee2a4a190d415a3664d893ad Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 1 May 2023 13:08:01 -0600 Subject: [PATCH 11/41] Small refactor --- src/KOKKOS/fix_langevin_kokkos.cpp | 9 ++++++ src/KOKKOS/fix_langevin_kokkos.h | 1 + src/KOKKOS/fix_nve_kokkos.cpp | 52 ++++++++++++++++++++++++------ src/KOKKOS/fix_nve_kokkos.h | 17 +++++----- src/KOKKOS/modify_kokkos.cpp | 4 +-- src/KOKKOS/modify_kokkos.h | 2 +- src/KOKKOS/verlet_kokkos.cpp | 2 +- src/fix.h | 2 +- src/modify.h | 2 +- 9 files changed, 67 insertions(+), 24 deletions(-) diff --git a/src/KOKKOS/fix_langevin_kokkos.cpp b/src/KOKKOS/fix_langevin_kokkos.cpp index b7305644c9..14eb0f1ab7 100644 --- a/src/KOKKOS/fix_langevin_kokkos.cpp +++ b/src/KOKKOS/fix_langevin_kokkos.cpp @@ -44,6 +44,7 @@ FixLangevinKokkos::FixLangevinKokkos(LAMMPS *lmp, int narg, char **a FixLangevin(lmp, narg, arg),rand_pool(seed + comm->me) { kokkosable = 1; + fuse_integrate_flag = 1; atomKK = (AtomKokkos *) atom; int ntypes = atomKK->ntypes; @@ -169,6 +170,14 @@ void FixLangevinKokkos::initial_integrate_item(int i) const /* ---------------------------------------------------------------------- */ +template +void FixLangevinKokkos::fused_integrate(int vflag) +{ + initial_integrate(vflag); +} + +/* ---------------------------------------------------------------------- */ + template void FixLangevinKokkos::post_force(int /*vflag*/) { diff --git a/src/KOKKOS/fix_langevin_kokkos.h b/src/KOKKOS/fix_langevin_kokkos.h index f7142e6286..0bd628270e 100644 --- a/src/KOKKOS/fix_langevin_kokkos.h +++ b/src/KOKKOS/fix_langevin_kokkos.h @@ -69,6 +69,7 @@ namespace LAMMPS_NS { void cleanup_copy(); void init() override; void initial_integrate(int) override; + void fused_integrate(int) override; void post_force(int) override; void reset_dt() override; void grow_arrays(int) override; diff --git a/src/KOKKOS/fix_nve_kokkos.cpp b/src/KOKKOS/fix_nve_kokkos.cpp index c26e26a02c..b8236c2657 100644 --- a/src/KOKKOS/fix_nve_kokkos.cpp +++ b/src/KOKKOS/fix_nve_kokkos.cpp @@ -160,21 +160,12 @@ void FixNVEKokkos::final_integrate_rmass_item(int i) const } } -/* ---------------------------------------------------------------------- */ - -template -void FixNVEKokkos::cleanup_copy() -{ - id = style = nullptr; - vatom = nullptr; -} - /* ---------------------------------------------------------------------- allow for both per-type and per-atom mass ------------------------------------------------------------------------- */ template -void FixNVEKokkos::fused_integrate() +void FixNVEKokkos::fused_integrate(int /*vflag*/) { atomKK->sync(execution_space,datamask_read); @@ -199,6 +190,47 @@ void FixNVEKokkos::fused_integrate() atomKK->modified(execution_space,datamask_modify); } +template +KOKKOS_INLINE_FUNCTION +void FixNVEKokkos::fused_integrate_item(int i) const +{ + if (mask[i] & groupbit) { + const double dtfm = 2.0 * dtf / mass[type[i]]; + v(i,0) += dtfm * f(i,0); + v(i,1) += dtfm * f(i,1); + v(i,2) += dtfm * f(i,2); + x(i,0) += dtv * v(i,0); + x(i,1) += dtv * v(i,1); + x(i,2) += dtv * v(i,2); + } +} + +template +KOKKOS_INLINE_FUNCTION +void FixNVEKokkos::fused_integrate_rmass_item(int i) const +{ + if (mask[i] & groupbit) { + const double dtfm = 2.0 * dtf / rmass[i]; + v(i,0) += dtfm * f(i,0); + v(i,1) += dtfm * f(i,1); + v(i,2) += dtfm * f(i,2); + x(i,0) += dtv * v(i,0); + x(i,1) += dtv * v(i,1); + x(i,2) += dtv * v(i,2); + } +} + +/* ---------------------------------------------------------------------- */ + +template +void FixNVEKokkos::cleanup_copy() +{ + id = style = nullptr; + vatom = nullptr; +} + +/* ---------------------------------------------------------------------- */ + namespace LAMMPS_NS { template class FixNVEKokkos; #ifdef LMP_KOKKOS_GPU diff --git a/src/KOKKOS/fix_nve_kokkos.h b/src/KOKKOS/fix_nve_kokkos.h index a10f846cba..a1e8e1398c 100644 --- a/src/KOKKOS/fix_nve_kokkos.h +++ b/src/KOKKOS/fix_nve_kokkos.h @@ -46,7 +46,7 @@ class FixNVEKokkos : public FixNVE { void init() override; void initial_integrate(int) override; void final_integrate() override; - void fused_integrate() override; + void fused_integrate(int) override; KOKKOS_INLINE_FUNCTION void initial_integrate_item(int) const; @@ -56,6 +56,10 @@ class FixNVEKokkos : public FixNVE { void final_integrate_item(int) const; KOKKOS_INLINE_FUNCTION void final_integrate_rmass_item(int) const; + KOKKOS_INLINE_FUNCTION + void fused_integrate_item(int) const; + KOKKOS_INLINE_FUNCTION + void fused_integrate_rmass_item(int) const; private: @@ -106,13 +110,10 @@ struct FixNVEKokkosFusedIntegrateFunctor { c(*c_ptr) {c.cleanup_copy();}; KOKKOS_INLINE_FUNCTION void operator()(const int i) const { - if (RMass) { - c.final_integrate_rmass_item(i); - c.initial_integrate_rmass_item(i); - } else { - c.final_integrate_item(i); - c.initial_integrate_item(i); - } + if (RMass) + c.fused_integrate_rmass_item(i); + else + c.fused_integrate_item(i); } }; diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp index 44fcb48727..0b81a1cabb 100644 --- a/src/KOKKOS/modify_kokkos.cpp +++ b/src/KOKKOS/modify_kokkos.cpp @@ -396,14 +396,14 @@ void ModifyKokkos::final_integrate() fused initial and final integrate call, only for relevant fixes ------------------------------------------------------------------------- */ -void ModifyKokkos::fused_integrate() +void ModifyKokkos::fused_integrate(int vflag) { for (int i = 0; i < n_final_integrate; i++) { atomKK->sync(fix[list_final_integrate[i]]->execution_space, fix[list_final_integrate[i]]->datamask_read); int prev_auto_sync = lmp->kokkos->auto_sync; if (!fix[list_final_integrate[i]]->kokkosable) lmp->kokkos->auto_sync = 1; - fix[list_final_integrate[i]]->fused_integrate(); + fix[list_final_integrate[i]]->fused_integrate(vflag); lmp->kokkos->auto_sync = prev_auto_sync; atomKK->modified(fix[list_final_integrate[i]]->execution_space, fix[list_final_integrate[i]]->datamask_modify); diff --git a/src/KOKKOS/modify_kokkos.h b/src/KOKKOS/modify_kokkos.h index 8dd9e6d9df..527518219c 100644 --- a/src/KOKKOS/modify_kokkos.h +++ b/src/KOKKOS/modify_kokkos.h @@ -39,7 +39,7 @@ class ModifyKokkos : public Modify { void pre_reverse(int,int) override; void post_force(int) override; void final_integrate() override; - void fused_integrate() override; + void fused_integrate(int) override; void end_of_step() override; double energy_couple() override; double energy_global() override; diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index bd15859644..3b88b34b38 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -498,7 +498,7 @@ void VerletKokkos::run(int n) if (n_post_force) modify->post_force(vflag); - if (fuse_integrate) modify->fused_integrate(); + if (fuse_integrate) modify->fused_integrate(vflag); else modify->final_integrate(); if (n_end_of_step) modify->end_of_step(); diff --git a/src/fix.h b/src/fix.h index 23b120d989..9676651afb 100644 --- a/src/fix.h +++ b/src/fix.h @@ -153,7 +153,6 @@ class Fix : protected Pointers { virtual void setup_pre_reverse(int, int) {} virtual void min_setup(int) {} virtual void initial_integrate(int) {} - virtual void fused_integrate() {} virtual void post_integrate() {} virtual void pre_exchange() {} virtual void pre_neighbor() {} @@ -162,6 +161,7 @@ class Fix : protected Pointers { virtual void pre_reverse(int, int) {} virtual void post_force(int) {} virtual void final_integrate() {} + virtual void fused_integrate(int) {} virtual void end_of_step() {} virtual void post_run() {} virtual void write_restart(FILE *) {} diff --git a/src/modify.h b/src/modify.h index 9686a9a7ec..6ca4b4ad26 100644 --- a/src/modify.h +++ b/src/modify.h @@ -61,7 +61,6 @@ class Modify : protected Pointers { virtual void setup_pre_force(int); virtual void setup_pre_reverse(int, int); virtual void initial_integrate(int); - virtual void fused_integrate() {} virtual void post_integrate(); virtual void pre_exchange(); virtual void pre_neighbor(); @@ -70,6 +69,7 @@ class Modify : protected Pointers { virtual void pre_reverse(int, int); virtual void post_force(int); virtual void final_integrate(); + virtual void fused_integrate(int) {} virtual void end_of_step(); virtual double energy_couple(); virtual double energy_global(); From 757130f9d049561b0a756c0b0b41216afbfbb572 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 1 May 2023 13:09:46 -0600 Subject: [PATCH 12/41] whitespace --- src/KOKKOS/fix_nve_kokkos.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KOKKOS/fix_nve_kokkos.cpp b/src/KOKKOS/fix_nve_kokkos.cpp index b8236c2657..fba8911f31 100644 --- a/src/KOKKOS/fix_nve_kokkos.cpp +++ b/src/KOKKOS/fix_nve_kokkos.cpp @@ -224,7 +224,7 @@ void FixNVEKokkos::fused_integrate_rmass_item(int i) const template void FixNVEKokkos::cleanup_copy() -{ +{ id = style = nullptr; vatom = nullptr; } From cd972884ff979e83daf1a8af494070db99246ced Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 1 May 2023 13:13:23 -0600 Subject: [PATCH 13/41] Clarify with comment --- src/KOKKOS/verlet_kokkos.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index 3b88b34b38..becdf67c9c 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -358,12 +358,15 @@ void VerletKokkos::run(int n) } } + // check if kernels can be fused, must come after initial_integrate + + fuse_check(i,n); + // force computations // important for pair to come before bonded contributions // since some bonded potentials tally pairwise energy/virial // and Pair:ev_tally() needs to be called before any tallying - fuse_check(i,n); if (!fuse_force_clear) force_clear(); From 76f9336de41939f459e5a597fdf9756a9574c0e9 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 1 May 2023 13:16:16 -0600 Subject: [PATCH 14/41] whitespace --- src/KOKKOS/fix_nve_kokkos.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/fix_nve_kokkos.cpp b/src/KOKKOS/fix_nve_kokkos.cpp index fba8911f31..59cc90c088 100644 --- a/src/KOKKOS/fix_nve_kokkos.cpp +++ b/src/KOKKOS/fix_nve_kokkos.cpp @@ -193,8 +193,8 @@ void FixNVEKokkos::fused_integrate(int /*vflag*/) template KOKKOS_INLINE_FUNCTION void FixNVEKokkos::fused_integrate_item(int i) const -{ - if (mask[i] & groupbit) { +{ + if (mask[i] & groupbit) { const double dtfm = 2.0 * dtf / mass[type[i]]; v(i,0) += dtfm * f(i,0); v(i,1) += dtfm * f(i,1); From f1b29941fdf1aa7a6da41c4899d07e933640e60d Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 1 May 2023 13:34:58 -0600 Subject: [PATCH 15/41] Port changes to nve/sphere --- src/KOKKOS/fix_nve_sphere_kokkos.cpp | 67 ++++++++++++++++++++++++++++ src/KOKKOS/fix_nve_sphere_kokkos.h | 14 ++++++ 2 files changed, 81 insertions(+) diff --git a/src/KOKKOS/fix_nve_sphere_kokkos.cpp b/src/KOKKOS/fix_nve_sphere_kokkos.cpp index 51e57b839e..ecea4f31d1 100644 --- a/src/KOKKOS/fix_nve_sphere_kokkos.cpp +++ b/src/KOKKOS/fix_nve_sphere_kokkos.cpp @@ -28,6 +28,7 @@ FixNVESphereKokkos::FixNVESphereKokkos(LAMMPS *lmp, int narg, char * FixNVESphere(lmp, narg, arg) { kokkosable = 1; + fuse_integrate_flag = 1; atomKK = (AtomKokkos *)atom; execution_space = ExecutionSpaceFromDevice::space; @@ -164,6 +165,72 @@ void FixNVESphereKokkos::final_integrate_item(const int i) const } } +/* ---------------------------------------------------------------------- */ + +template +void FixNVESphereKokkos::fused_integrate(int /*vflag*/) +{ + if (extra == DIPOLE) + atomKK->sync(execution_space, X_MASK | V_MASK | OMEGA_MASK| F_MASK | TORQUE_MASK | RMASS_MASK | RADIUS_MASK | MASK_MASK | MU_MASK); + else + atomKK->sync(execution_space, X_MASK | V_MASK | OMEGA_MASK| F_MASK | TORQUE_MASK | RMASS_MASK | RADIUS_MASK | MASK_MASK); + + v = atomKK->k_v.view(); + omega = atomKK->k_omega.view(); + f = atomKK->k_f.view(); + torque = atomKK->k_torque.view(); + mask = atomKK->k_mask.view(); + rmass = atomKK->k_rmass.view(); + radius = atomKK->k_radius.view(); + mu = atomKK->k_mu.view(); + + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + FixNVESphereKokkosFusedIntegrateFunctor f(this); + Kokkos::parallel_for(nlocal,f); + + if (extra == DIPOLE) + atomKK->modified(execution_space, X_MASK | V_MASK | OMEGA_MASK | MU_MASK); + else + atomKK->modified(execution_space, X_MASK | V_MASK | OMEGA_MASK); +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void FixNVESphereKokkos::fused_integrate_item(const int i) const +{ + const double dtfrotate = dtf / inertia; + + if (mask(i) & groupbit) { + const double dtfm = 2.0 * dtf / rmass(i); + v(i,0) += dtfm * f(i,0); + v(i,1) += dtfm * f(i,1); + v(i,2) += dtfm * f(i,2); + x(i,0) += dtv * v(i,0); + x(i,1) += dtv * v(i,1); + x(i,2) += dtv * v(i,2); + + const double dtirotate = 2.0 * dtfrotate / (radius(i)*radius(i)*rmass(i)); + omega(i,0) += dtirotate * torque(i,0); + omega(i,1) += dtirotate * torque(i,1); + omega(i,2) += dtirotate * torque(i,2); + + if (extra == DIPOLE) { + const double g0 = mu(i,0) + dtv * (omega(i,1) * mu(i,2) - omega(i,2) * mu(i,1)); + const double g1 = mu(i,1) + dtv * (omega(i,2) * mu(i,0) - omega(i,0) * mu(i,2)); + const double g2 = mu(i,2) + dtv * (omega(i,0) * mu(i,1) - omega(i,1) * mu(i,0)); + const double msq = g0*g0 + g1*g1 + g2*g2; + const double scale = mu(i,3)/sqrt(msq); + mu(i,0) = g0*scale; + mu(i,1) = g1*scale; + mu(i,2) = g2*scale; + } + } +} + namespace LAMMPS_NS { template class FixNVESphereKokkos; #ifdef LMP_KOKKOS_GPU diff --git a/src/KOKKOS/fix_nve_sphere_kokkos.h b/src/KOKKOS/fix_nve_sphere_kokkos.h index 268b4ea9ce..02acb466f2 100644 --- a/src/KOKKOS/fix_nve_sphere_kokkos.h +++ b/src/KOKKOS/fix_nve_sphere_kokkos.h @@ -37,11 +37,14 @@ class FixNVESphereKokkos : public FixNVESphere { void init() override; void initial_integrate(int) override; void final_integrate() override; + void fused_integrate(int) override; KOKKOS_INLINE_FUNCTION void initial_integrate_item(const int i) const; KOKKOS_INLINE_FUNCTION void final_integrate_item(const int i) const; + KOKKOS_INLINE_FUNCTION + void fused_integrate_item(int) const; private: typename ArrayTypes::t_x_array x; @@ -77,6 +80,17 @@ struct FixNVESphereKokkosFinalIntegrateFunctor { } }; +template +struct FixNVESphereKokkosFusedIntegrateFunctor { + typedef DeviceType device_type; + FixNVESphereKokkos c; + FixNVESphereKokkosFusedIntegrateFunctor(FixNVESphereKokkos *c_ptr): c(*c_ptr) { c.cleanup_copy(); } + KOKKOS_INLINE_FUNCTION + void operator()(const int i) const { + c.fused_integrate_item(i); + } +}; + } // namespace LAMMPS_NS #endif // LMP_FIX_NVE_SPHERE_KOKKOS_H From ef9ce62aa8f89b9ce067cda9a1b8fc6f6bb78cd4 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 1 May 2023 15:17:19 -0600 Subject: [PATCH 16/41] update broken bond doc page --- doc/src/Howto_broken_bonds.rst | 75 ++++++++++++++++++---------------- 1 file changed, 40 insertions(+), 35 deletions(-) diff --git a/doc/src/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index 88c8acb87f..4fd0bcc30f 100644 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -1,48 +1,53 @@ Broken Bonds ============ -Typically, bond interactions persist for the duration of a simulation in -LAMMPS. However, there are some exceptions that allow for bonds to -break, including the :doc:`quartic bond style ` and the -bond styles in the :doc:`BPM package ` which contains the -:doc:`bpm/spring ` and :doc:`bpm/rotational -` bond styles. In these cases, a bond can be broken -if it is stretched beyond a user-defined threshold. LAMMPS accomplishes -this by setting the bond type to 0, such that the bond force is no -longer computed. +Typically, bond interactions persist for the duration of a simulation +in LAMMPS. However, some commands allow bonds to break, including the +following: -Users are normally able to weight the contribution of pair forces to atoms -that are bonded using the :doc:`special_bonds command `. -When bonds break, this is not always the case. For the quartic bond style, -pair forces are always turned off between bonded particles. LAMMPS does -this via a computational sleight-of-hand. It subtracts the pairwise -interaction as part of the bond computation. When the bond breaks, the -subtraction stops. For this to work, the pairwise interaction must always -be computed by the :doc:`pair_style ` command, whether the bond -is broken or not. This means that :doc:`special_bonds ` must -be set to 1,1,1. After the bond breaks, the pairwise interaction between the -two atoms is turned on, since they are no longer bonded. +* :doc:`bond_style quartic ` +* :doc:`fix bond/break ` +* :doc:`fix bond/react ` +* :doc:`BPM package ` bond styles -In the BPM package, one can either turn off all pair interactions between -bonded particles or leave them on, overlaying pair forces on top of bond -forces. To remove pair forces, the special bond list is dynamically -updated. More details can be found on the :doc:`Howto BPM ` -page. +A bond can break if it is stretched beyond a user-defined threshold or +more generally if other criteria are met. -Bonds can also be broken by fixes which change bond topology, including -:doc:`fix bond/break ` and -:doc:`fix bond/react `. These fixes will automatically -trigger a rebuild of the neighbor list and update special bond data structures -when bonds are broken. +For the quartic bond style, when a bond is broken its bond type is set +to 0 and pairwise forces between the two atoms in the broken bond are +"turned on". Angles, dihedrals, etc cannot be defined for the system +when :doc:`bond_style quartic ` is used. -Note that when bonds are dumped to a file via the :doc:`dump local ` command, bonds with type 0 are not included. The -:doc:`delete_bonds ` command can also be used to query the -status of broken bonds or permanently delete them, e.g.: +The :doc:`fix bond/break ` and :doc:`fix bond/react +` commands allow breaking of bonds within a molecular +topology with also defines angles, dihedrals, etc. These fixes will +update internal topology data structures when bonds are broken, so +that the appropriate angle, dihederal, etc interactions are also +turned off. They will also trigger a rebuild of the neighbor list +when this occurs, to turn on the appropriate pairwise forces. + +In the BPM package, one can either turn off all pair interactions +between bonded particles or leave them on, overlaying pair forces on +top of bond forces. To remove pair forces, the special bond list is +dynamically updated. More details can be found on the :doc:`Howto BPM +` page. + +Note that when bonds are dumped to a file via the :doc:`dump local +` command, bonds with type 0 are not included. + +The :doc:`delete_bonds ` command can also be used to +query the status of broken bonds (type = 0) or permanently delete +them, e.g.: .. code-block:: LAMMPS delete_bonds all stats delete_bonds all bond 0 remove -The compute :doc:`nbond/atom ` can also be used -to tally the current number of bonds per atom, excluding broken bonds. +The compute :doc:`count/type bond ` command +tallies the current number of bonds for each bond type. It also +tallies broken bonds with type = 0. + +The compute :doc:`nbond/atom ` command tallies the +current number of bonds each atom is part of, excluding broken bonds +with type = 0. From fae0ef6cdea8a531bb9a5c77bb8e32c8ab8284da Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 1 May 2023 15:21:56 -0600 Subject: [PATCH 17/41] update count/type doc page --- doc/src/compute_count_type.rst | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index ca56648297..9ab0b1009f 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -25,16 +25,22 @@ Examples Description """"""""""" -Define a computation that counts the current number of atoms by atom -type or the number of bonds by bond type. The latter can be useful in -reactive simulations where bonds are broken or created. +Define a computation that counts the current number of atoms for each +atom type or the number of bonds for each bond type. The former can +be useful if atoms are added to or deleted from the system in random +ways, e.g. via the :doc:`fix deposit `, :doc:`fix pour +`, or :doc:`fix evaporate ` commands. The +latter can be useful in reactive simulations where molecular bonds are +broken or created. -Note that for this command, bonds are the topological ones enumerated -in a data file, initially read by the :doc:`read_data ` -command. They do not refer to bonds defined on-the-fly by bond-order -or reactive pair styles. +Note that for this command, bonds are defined as the topological kind +enumerated in a data file, initially read by the :doc:`read_data +` command or defined by the :doc:`molecule ` +command. They do not refer to implicit bonds defined on-the-fly by +bond-order or reactive pair styles based on the current conformation +of small clusters of atoms. -These commands can create and break toplogical bonds: +These commands can create and/or break topological bonds: * :doc:`fix bond/react ` * :doc:`fix bond/create ` @@ -50,9 +56,11 @@ type is tallied. Only bonds with both atoms in the specified group are counted. For {mode} = {bond}, broken bonds with a bond type of zero are also -counted. See the :doc:`Howto broken bonds ` doc -page for details. Note that the group setting is ignored for broken -bonds; all broken bonds in the system are counted. +counted. Some of the commands listed above which break bonds, do this +by setting their types to zero. See the :doc:`Howto broken bonds +` doc page for details. Note that the group +setting is ignored for broken bonds; all broken bonds in the system +are counted. ---------- From f8058f039fe77e30974d2fb172cebed33e68f38b Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 1 May 2023 15:26:07 -0600 Subject: [PATCH 18/41] Fix bug --- src/KOKKOS/fix_nve_sphere_kokkos.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/KOKKOS/fix_nve_sphere_kokkos.cpp b/src/KOKKOS/fix_nve_sphere_kokkos.cpp index ecea4f31d1..38f6a40792 100644 --- a/src/KOKKOS/fix_nve_sphere_kokkos.cpp +++ b/src/KOKKOS/fix_nve_sphere_kokkos.cpp @@ -175,6 +175,7 @@ void FixNVESphereKokkos::fused_integrate(int /*vflag*/) else atomKK->sync(execution_space, X_MASK | V_MASK | OMEGA_MASK| F_MASK | TORQUE_MASK | RMASS_MASK | RADIUS_MASK | MASK_MASK); + x = atomKK->k_x.view(); v = atomKK->k_v.view(); omega = atomKK->k_omega.view(); f = atomKK->k_f.view(); From dbcc19aef4ab0386f6e1ba7e011926ca2e190cc7 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 1 May 2023 15:42:17 -0600 Subject: [PATCH 19/41] Fix uninitialized variables --- src/KOKKOS/verlet_kokkos.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index becdf67c9c..7570f1d8fa 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -276,6 +276,9 @@ void VerletKokkos::run(int n) lmp->kokkos->auto_sync = 0; + fuse_integrate = 0; + fuse_force_clear = 0; + if (atomKK->sortfreq > 0) sortflag = 1; else sortflag = 0; From 478058119b051ea256a519f843abecc8c9a7bb43 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 29 Apr 2023 02:59:04 -0400 Subject: [PATCH 20/41] integrate CMake build procedure for tools/phonon --- cmake/CMakeLists.txt | 2 +- cmake/Modules/Tools.cmake | 2 + tools/phonon/CMakeLists.spglib | 54 +++++ tools/phonon/CMakeLists.txt | 99 ++++++++ tools/phonon/Makefile | 67 ------ tools/phonon/README | 56 +++-- tools/phonon/disp.cpp | 43 ++-- tools/phonon/dynmat.cpp | 72 ++++-- tools/phonon/dynmat.h | 20 +- tools/phonon/global.h | 2 +- tools/phonon/green.cpp | 13 +- tools/phonon/green.h | 5 +- tools/phonon/input.cpp | 35 +++ tools/phonon/input.h | 17 ++ tools/phonon/interpolate.cpp | 21 +- tools/phonon/interpolate.h | 15 +- tools/phonon/kpath.cpp | 39 ++- tools/phonon/kpath.h | 17 +- tools/phonon/main.cpp | 4 - tools/phonon/memory.cpp | 6 +- tools/phonon/memory.h | 308 +++++------------------- tools/phonon/phonon.cpp | 176 +++++++++----- tools/phonon/phonon.h | 16 +- tools/phonon/phonopy.cpp | 141 ++++++----- tools/phonon/phonopy.h | 16 +- tools/phonon/timer.cpp | 3 + tools/phonon/timer.h | 4 +- tools/phonon/tricubic/CMakeLists.txt | 18 ++ tools/phonon/tricubic/LICENSE | 340 +++++++++++++++++++++++++++ tools/phonon/tricubic/README.md | 28 +++ tools/phonon/tricubic/tricubic.cpp | 185 +++++++++++++++ tools/phonon/tricubic/tricubic.h | 10 + tools/phonon/version.h | 1 - tools/phonon/zheevd.h | 16 ++ 34 files changed, 1260 insertions(+), 591 deletions(-) create mode 100644 tools/phonon/CMakeLists.spglib create mode 100644 tools/phonon/CMakeLists.txt delete mode 100644 tools/phonon/Makefile create mode 100644 tools/phonon/input.cpp create mode 100644 tools/phonon/input.h create mode 100644 tools/phonon/tricubic/CMakeLists.txt create mode 100644 tools/phonon/tricubic/LICENSE create mode 100644 tools/phonon/tricubic/README.md create mode 100644 tools/phonon/tricubic/tricubic.cpp create mode 100644 tools/phonon/tricubic/tricubic.h delete mode 100644 tools/phonon/version.h create mode 100644 tools/phonon/zheevd.h diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index a71347c2c4..88b7ec422e 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -440,7 +440,7 @@ if(BUILD_OMP) target_link_libraries(lmp PRIVATE OpenMP::OpenMP_CXX) endif() -if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE) +if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE OR BUILD_TOOLS) enable_language(C) if (NOT USE_INTERNAL_LINALG) find_package(LAPACK) diff --git a/cmake/Modules/Tools.cmake b/cmake/Modules/Tools.cmake index c4c33cc40c..ba91a557bc 100644 --- a/cmake/Modules/Tools.cmake +++ b/cmake/Modules/Tools.cmake @@ -63,6 +63,8 @@ if(BUILD_LAMMPS_SHELL) install(TARGETS lammps-shell EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR}) install(DIRECTORY ${LAMMPS_TOOLS_DIR}/lammps-shell/icons DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/) install(FILES ${LAMMPS_TOOLS_DIR}/lammps-shell/lammps-shell.desktop DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/applications/) + + add_subdirectory(${LAMMPS_TOOLS_DIR}/phonon ${CMAKE_BINARY_DIR}/phana_build) endif() diff --git a/tools/phonon/CMakeLists.spglib b/tools/phonon/CMakeLists.spglib new file mode 100644 index 0000000000..9e22445b80 --- /dev/null +++ b/tools/phonon/CMakeLists.spglib @@ -0,0 +1,54 @@ +cmake_minimum_required(VERSION 3.10) +project(spglib C) +set(CMAKE_MACOSX_RPATH 1) +set(CMAKE_C_FLAGS_RELEASE "-Wall -O2") +set(CMAKE_C_FLAGS_DEBUG "-g -DSPGDEBUG -DSPGWARNING") +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif(NOT CMAKE_BUILD_TYPE) +message(STATUS "Build type: ${CMAKE_BUILD_TYPE}") + +set(CMAKE_POSITION_INDEPENDENT_CODE ON) + + +# Version numbers +file(READ ${PROJECT_SOURCE_DIR}/src/version.h version_file) +string(REGEX MATCH "SPGLIB_MAJOR_VERSION ([0-9]+)" spglib_major_version ${version_file}) +set(spglib_major_version ${CMAKE_MATCH_1}) +string(REGEX MATCH "SPGLIB_MINOR_VERSION ([0-9]+)" spglib_minor_version ${version_file}) +set(spglib_minor_version ${CMAKE_MATCH_1}) +string(REGEX MATCH "SPGLIB_MICRO_VERSION ([0-9]+)" spglib_micro_version ${version_file}) +set(spglib_micro_version ${CMAKE_MATCH_1}) +set(serial "${spglib_major_version}.${spglib_minor_version}.${spglib_micro_version}") +set(soserial "1") + +# Source code +include_directories("${PROJECT_SOURCE_DIR}/src") +set(SOURCES ${PROJECT_SOURCE_DIR}/src/arithmetic.c + ${PROJECT_SOURCE_DIR}/src/cell.c + ${PROJECT_SOURCE_DIR}/src/debug.c + ${PROJECT_SOURCE_DIR}/src/delaunay.c + ${PROJECT_SOURCE_DIR}/src/determination.c + ${PROJECT_SOURCE_DIR}/src/hall_symbol.c + ${PROJECT_SOURCE_DIR}/src/kgrid.c + ${PROJECT_SOURCE_DIR}/src/kpoint.c + ${PROJECT_SOURCE_DIR}/src/mathfunc.c + ${PROJECT_SOURCE_DIR}/src/niggli.c + ${PROJECT_SOURCE_DIR}/src/overlap.c + ${PROJECT_SOURCE_DIR}/src/pointgroup.c + ${PROJECT_SOURCE_DIR}/src/primitive.c + ${PROJECT_SOURCE_DIR}/src/refinement.c + ${PROJECT_SOURCE_DIR}/src/site_symmetry.c + ${PROJECT_SOURCE_DIR}/src/sitesym_database.c + ${PROJECT_SOURCE_DIR}/src/spacegroup.c + ${PROJECT_SOURCE_DIR}/src/spg_database.c + ${PROJECT_SOURCE_DIR}/src/spglib.c + ${PROJECT_SOURCE_DIR}/src/spin.c + ${PROJECT_SOURCE_DIR}/src/symmetry.c) + +add_library(symspg STATIC ${SOURCES}) +install(TARGETS symspg LIBRARY) + +# Header file +install(FILES ${PROJECT_SOURCE_DIR}/src/spglib.h DESTINATION ${CMAKE_INSTALL_PREFIX}/include) + diff --git a/tools/phonon/CMakeLists.txt b/tools/phonon/CMakeLists.txt new file mode 100644 index 0000000000..283b995590 --- /dev/null +++ b/tools/phonon/CMakeLists.txt @@ -0,0 +1,99 @@ + +# Support Linux from Ubuntu 20.04LTS onward, CentOS 7.x (with EPEL), +# macOS, MSVC 2019 (=Version 16) +cmake_minimum_required(VERSION 3.16) + +# set timestamp of downloaded files to that of archive +if(POLICY CMP0135) + cmake_policy(SET CMP0135 NEW) +endif() + +# set up project +set(PHANA_MINOR_VERSION 48) +project(phonon VERSION ${PHANA_MINOR_VERSION} + DESCRIPTION "Fix phonon post-processor" + LANGUAGES CXX C) +set(CMAKE_POSITION_INDEPENDENT_CODE ON) +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE RelWithDebInfo) +endif() + +# hacks for MSVC to prevent lots of pointless warnings about "unsafe" functions, +# padding and Spectre mitigation +if(MSVC) + add_compile_options(/wd4244) + add_compile_options(/wd4267) + add_compile_options(/wd4711) + add_compile_options(/wd4820) + add_compile_options(/wd5045) + add_compile_definitions(_CRT_SECURE_NO_WARNINGS) +endif() + +set(CMAKE_PROJECT_VERSION ${PHANA_MINOR_VERSION}) +configure_file(version.h.in version.h @ONLY) +add_executable(phana + main.cpp + disp.cpp + dynmat.cpp + green.cpp + input.cpp + interpolate.cpp + kpath.cpp + memory.cpp + phonon.cpp + phonopy.cpp + qnodes.cpp + timer.cpp +) +target_include_directories(phana PUBLIC $) + +find_package(FFTW3) +if(FFTW3_FOUND) + target_compile_definitions(phana PRIVATE FFTW3) + target_link_libraries(phana PRIVATE FFTW3::FFTW3) +endif() + +# build bundeled libraries +add_subdirectory(tricubic) + +option(USE_SPGLIB "Download and use spglib for phonon DOS and other optional properties" ON) +if(USE_SPGLIB) + set(SPGLIB_URL "https://github.com/spglib/spglib/archive/refs/tags/v1.11.2.1.tar.gz" CACHE STRING "URL for spglib v1.x tarball") + set(SPGLIB_MD5 "3089782bc85b5034dd4765a18ee70bc7" CACHE STRING "MD5 checksum for spglib tarball") + mark_as_advanced(SPGLIB_URL) + mark_as_advanced(SPGLIB_MD5) + GetFallbackURL(SPGLIB_URL SPGLIB_FALLBACK) + + include(ExternalProject) + ExternalProject_Add(spglib_build + URL ${SPGLIB_URL} ${SPGLIB_FALLBACK} + URL_MD5 ${SPGLIB_MD5} + PREFIX ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext + CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext + -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} + -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} + -DCMAKE_MAKE_PROGRAM=${CMAKE_MAKE_PROGRAM} + -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE} + -DCMAKE_POSITION_INDEPENDENT_CODE=ON + UPDATE_COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CMAKE_CURRENT_SOURCE_DIR}/CMakeLists.spglib ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/src/spglib_build/CMakeLists.txt + INSTALL_COMMAND ${CMAKE_COMMAND} --build ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/src/spglib_build-build --target install + BUILD_BYPRODUCTS "${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib/${CMAKE_STATIC_LIBRARY_PREFIX}symspg${CMAKE_STATIC_LIBRARY_SUFFIX}" + ) + + # workaround for older CMake versions + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib) + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/include) + + add_library(SPGLIB::SYMSPG UNKNOWN IMPORTED) + add_dependencies(SPGLIB::SYMSPG spglib_build) + set_target_properties(SPGLIB::SYMSPG PROPERTIES + IMPORTED_LOCATION "${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib/${CMAKE_STATIC_LIBRARY_PREFIX}symspg${CMAKE_STATIC_LIBRARY_SUFFIX}" + INTERFACE_INCLUDE_DIRECTORIES ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/include + ) + + target_compile_definitions(phana PRIVATE UseSPG) + target_link_libraries(phana PRIVATE SPGLIB::SYMSPG) +endif() + +target_link_libraries(phana PRIVATE tricubic ${LAPACK_LIBRARIES}) +install(TARGETS phana EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/tools/phonon/Makefile b/tools/phonon/Makefile deleted file mode 100644 index 6afae5f312..0000000000 --- a/tools/phonon/Makefile +++ /dev/null @@ -1,67 +0,0 @@ -.SUFFIXES : .o .cpp -# compiler and flags -CC = g++ -Wno-unused-result -LINK = $(CC) -static -CFLAGS = -O3 $(DEBUG) $(UFLAG) - -OFLAGS = -O3 $(DEBUG) -INC = $(LPKINC) $(TCINC) $(SPGINC) $(FFTINC) -LIB = $(LPKLIB) $(TCLIB) $(SPGLIB) $(FFTLIB) - -# cLapack library needed -LPKINC = -I/opt/clapack/3.2.1/include -LPKLIB = -L/opt/clapack/3.2.1/lib -lclapack -lblas -lf2c -lm - -# Tricubic library needed -TCINC = -I/opt/tricubic/1.0/include -TCLIB = -L/opt/tricubic/1.0/lib -ltricubic - -# spglib, used to get the irreducible q-points -# if SFLAG is not set, spglib won't be used. -SFLAG = -DUseSPG -SPGINC = -I/opt/spglib/1.9.7/include/spglib -SPGLIB = -L/opt/spglib/1.9.7/lib -lsymspg - -# FFTW 3， used to deduce the force constants in real space -# if FFLAG is not set, fftw won't be used. -FFLAG = -DFFTW3 -FFTINC = -I/opt/fftw/3.3.7/include -FFTLIB = -L/opt/fftw/3.3.7/lib -lfftw3 - -# Debug flags -# DEBUG = -g -DDEBUG -UFLAG = $(SFLAG) $(FFLAG) - -#==================================================================== -ROOT = phana -# executable name -EXE = $(ROOT) -#==================================================================== -# source and rules -SRC = $(wildcard *.cpp) -OBJ = $(SRC:.cpp=.o) - -#==================================================================== -all: ${EXE} - -${EXE}: $(OBJ) - $(LINK) $(OFLAGS) $(OBJ) $(LIB) -o $@ - -clean: - rm -f *.o *~ *.mod ${EXE} - -tar: - rm -f ${ROOT}.tar.gz; tar -czvf ${ROOT}.tar.gz *.cpp *.h Makefile README - -ver: - @echo "#define VERSION `git log|grep '^commit'|wc -l`" > version.h - -#==================================================================== -.f.o: - $(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $< -.f90.o: - $(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $< -.c.o: - $(CC) $(CFLAGS) -c $< -.cpp.o: - $(CC) $(CFLAGS) $(INC) -c $< diff --git a/tools/phonon/README b/tools/phonon/README index db69ac50c8..c1430f5940 100644 --- a/tools/phonon/README +++ b/tools/phonon/README @@ -5,34 +5,38 @@ analyse the phonon related information. #------------------------------------------------------------------------------- 1. Dependencies - The clapack library is needed to solve the eigen problems, - which could be downloaded from: - http://www.netlib.org/clapack/ - - The tricubic library is also needed to do tricubic interpolations, - which could now be obtained from: - https://github.com/nbigaouette/libtricubic/ - + The ZHEEVD LAPACK function is needed to solve the eigen problems. + A C++ compilable version based on CLAPACK is included in the linalg folder + and will be automatically built. + + The tricubic library is also needed to do tricubic interpolations. + A copy is included and will be automatically built. + The spglib is optionally needed, enabling one to evaluate the phonon density of states or vibrational thermal properties using only the irreducible q-points in the first Brillouin zone, as well as to evaluate the phonon dispersion curvers with the - automatic mode. Currently, the 1.8.3 version of spglib is used. - It can be obtained from: - http://spglib.sourceforge.net/ + automatic mode. Currently, version 1.11.2.1 of spglib is used. + It is automatically downloaded and compiled unless the -DUSE_SPGLIB=off + variable is set during CMake configuration. FFTW 3 might also be needed if you would like to interface with phonopy: necessary input files for phonopy will be prepared so that you can make use of the functions provided by phonopy. - FFTW 3 can be downloaded from: - http://www.fftw.org - + It is autodetected and used if available. + + FFTW 3 can be downloaded from: http://www.fftw.org + 2. Compilation - To compile the code, one needs therefore to install the above - libraries and set the paths correctly in the Makefile. - Once this is done, by typing - make - will yield the executable "phana". + To compile the code, one needs to have CMake version 3.16 + or later installed. + + The CMake configuration is done with: + cmake -S . -B build + And compilation then performed with: + cmake --build build + The phana (or phana.exe) executable is then available in + the "build" folder 3. Unit system The units of the output frequencies by this code is THz for @@ -46,6 +50,18 @@ 5. Bug report If any bug found, please drop a line to: konglt(at)sjtu.edu.cn +6. Precompiled executable + The "precompiled" folder contains a precompiled and statically + linked Linux executable for x86_64 CPUs. It should work on *any* + Linux machine with using the x86_64 architecture. It includes + spglib support but not fftw3. + +7. Portability + Build and use of phana has been successfully tested on: + - Fedora Linux 38 using GCC, Clang, and MinGW Linux2Windows cross-compiler + - macOS 12 (Monterey) using Xcode + - Windows 11 using Visual Studio 2022 with MSVC and Clang + #------------------------------------------------------------------------------- Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn -May 2020 +Aug 2021 diff --git a/tools/phonon/disp.cpp b/tools/phonon/disp.cpp index 8a53873383..79a22aeee4 100644 --- a/tools/phonon/disp.cpp +++ b/tools/phonon/disp.cpp @@ -1,10 +1,18 @@ -#include "string.h" -#include "qnodes.h" -#include "global.h" + #include "phonon.h" -#include "green.h" -#include "timer.h" + +#include "dynmat.h" +#include "global.h" +#include "input.h" #include "kpath.h" +#include "qnodes.h" + +#include +#include +#include +#include +#include +#include /*------------------------------------------------------------------------------ * Private method to evaluate the phonon dispersion curves @@ -13,19 +21,22 @@ void Phonon::pdisp() { // ask the output file name and write the header. char str[MAXLINE]; - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + puts("================================================================================"); + #ifdef UseSPG // ask method to generate q-lines int method = 2; printf("Please select your method to generate the phonon dispersion:\n"); printf(" 1. Manual, should always work;\n"); printf(" 2. Automatic, works only for 3D crystals (CMS49-299).\nYour choice [2]: "); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) method = atoi(strtok(str," \t\n\r\f")); method = 2 - method%2; printf("Your selection: %d\n", method); #endif printf("\nPlease input the filename to output the dispersion data [pdisp.dat]:"); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdisp.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "pdisp.dat"); char *ptr = strtok(str," \t\n\r\f"); char *fname = new char[strlen(ptr)+1]; strcpy(fname,ptr); @@ -45,9 +56,9 @@ void Phonon::pdisp() while (1){ for (int i = 0; i < 3; ++i) qstr[i] = qend[i]; - int quit = 0; printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]); - int n = count_words(fgets(str, MAXLINE, stdin)); + input->read_stdin(str); + int n = count_words(str); ptr = strtok(str, " \t\n\r\f"); if ((n == 1) && (strcmp(ptr,"q") == 0)) break; else if (n >= 3){ @@ -56,14 +67,18 @@ void Phonon::pdisp() qstr[2] = atof(strtok(NULL, " \t\n\r\f")); } - do printf("Please input the end q-point in unit of B1->B3: "); - while (count_words(fgets(str, MAXLINE, stdin)) < 3); + while ( 1 ){ + printf("Please input the end q-point in unit of B1->B3: "); + input->read_stdin(str); + if (count_words(str) >= 3) break; + } qend[0] = atof(strtok(str, " \t\n\r\f")); qend[1] = atof(strtok(NULL, " \t\n\r\f")); qend[2] = atof(strtok(NULL, " \t\n\r\f")); printf("Please input the # of points along the line [%d]: ", nq); - if (count_words(fgets(str, MAXLINE, stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) nq = atoi(strtok(str," \t\n\r\f")); nq = MAX(nq,2); double *qtmp = new double [3]; @@ -147,7 +162,7 @@ void Phonon::pdisp() printf("\nPhonon dispersion data are written to: %s, you can visualize the results\n", fname); printf("by invoking: `gnuplot pdisp.gnuplot; gv pdisp.eps`\n"); } - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + puts("================================================================================"); delete []fname; delete qnodes; diff --git a/tools/phonon/dynmat.cpp b/tools/phonon/dynmat.cpp index a6c4105547..aeb8fe6430 100644 --- a/tools/phonon/dynmat.cpp +++ b/tools/phonon/dynmat.cpp @@ -1,7 +1,16 @@ + #include "dynmat.h" -#include "math.h" -#include "version.h" + #include "global.h" +#include "input.h" +#include "interpolate.h" +#include "memory.h" +#include "version.h" +#include "zheevd.h" + +#include +#include +#include /* ---------------------------------------------------------------------------- * Class DynMat stores the Dynamic Matrix read from the binary file from @@ -9,6 +18,7 @@ * ---------------------------------------------------------------------------- */ DynMat::DynMat(int narg, char **arg) { + input = NULL; attyp = NULL; memory = NULL; M_inv_sqrt = NULL; @@ -19,6 +29,8 @@ DynMat::DynMat(int narg, char **arg) attyp = NULL; basis = NULL; flag_reset_gamma = flag_skip = 0; + symprec = -1.; + int flag_save = 0; // analyze the command line options int iarg = 1; @@ -29,9 +41,16 @@ DynMat::DynMat(int narg, char **arg) } else if (strcmp(arg[iarg], "-r") == 0){ flag_reset_gamma = 1; + } else if (strcmp(arg[iarg], "-p") == 0){ + if (++iarg >= narg) help(); + else symprec = fabs(atof(arg[iarg])); + } else if (strcmp(arg[iarg], "-h") == 0){ help(); + } else if (strcmp(arg[iarg], "-save") == 0){ + flag_save = 1; + } else { if (binfile) delete []binfile; int n = strlen(arg[iarg]) + 1; @@ -43,6 +62,8 @@ DynMat::DynMat(int narg, char **arg) } ShowVersion(); + input = new UserInput(flag_save); + // get the binary file name from user input if not found in command line char str[MAXLINE]; if (binfile == NULL) { @@ -50,7 +71,7 @@ DynMat::DynMat(int narg, char **arg) printf("\n"); do { printf("Please input the binary file name from fix_phonon: "); - fgets(str,MAXLINE,stdin); + input->read_stdin(str); ptr = strtok(str, " \n\t\r\f"); } while (ptr == NULL); @@ -137,17 +158,17 @@ DynMat::DynMat(int narg, char **arg) fclose(fp); exit(3); } - if (fread(basis[0], sizeof(double), fftdim, fp) != fftdim){ + if (fread(basis[0], sizeof(double), fftdim, fp) != (size_t)fftdim){ printf("\nError while reading basis info from file: %s\n", binfile); fclose(fp); exit(3); } - if (fread(&attyp[0], sizeof(int), nucell, fp) != nucell){ + if (fread(&attyp[0], sizeof(int), nucell, fp) != (size_t)nucell){ printf("\nError while reading atom types from file: %s\n", binfile); fclose(fp); exit(3); } - if (fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != nucell){ + if (fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != (size_t)nucell){ printf("\nError while reading atomic masses from file: %s\n", binfile); fclose(fp); exit(3); @@ -159,6 +180,7 @@ DynMat::DynMat(int narg, char **arg) // initialize interpolation interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all); + interpolate->input = input; if (flag_reset_gamma) interpolate->reset_gamma(); // Enforcing Austic Sum Rule @@ -217,7 +239,7 @@ void DynMat::writeDMq(double *q) printf("\n"); while ( 1 ){ printf("Please input the filename to output the DM at selected q: "); - fgets(str,MAXLINE,stdin); + input->read_stdin(str); ptr = strtok(str, " \r\t\n\f"); if (ptr) break; } @@ -264,9 +286,9 @@ void DynMat::writeDMq(double *q, const double qr, FILE *fp) int DynMat::geteigen(double *egv, int flag) { char jobz, uplo; - integer n, lda, lwork, lrwork, *iwork, liwork, info; + int n, lda, lwork, lrwork, *iwork, liwork, info; doublecomplex *work; - doublereal *w = &egv[0], *rwork; + double *w = &egv[0], *rwork; n = fftdim; if (flag) jobz = 'V'; @@ -348,7 +370,7 @@ void DynMat::EnforceASR() char str[MAXLINE]; int nasr = 20; if (nucell <= 1) nasr = 1; - printf("\n"); for (int i = 0; i < 80; ++i) printf("="); + printf("\n================================================================================"); // compute and display eigenvalues of Phi at gamma before ASR if (nucell > 100){ @@ -356,7 +378,7 @@ void DynMat::EnforceASR() fflush(stdout); } - double egvs[fftdim]; + double *egvs = new double[fftdim]; for (int i = 0; i < fftdim; ++i) for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j]; geteigen(egvs, 0); @@ -370,11 +392,11 @@ void DynMat::EnforceASR() // ask for iterations to enforce ASR printf("Please input the # of iterations to enforce ASR [%d]: ", nasr); - fgets(str,MAXLINE,stdin); + input->read_stdin(str); char *ptr = strtok(str," \t\n\r\f"); if (ptr) nasr = atoi(ptr); if (nasr < 1){ - for (int i=0; i<80; i++) printf("="); printf("\n"); + puts("================================================================================"); return; } @@ -439,9 +461,8 @@ void DynMat::EnforceASR() if (i%10 == 9) printf("\n"); if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;} } - printf("\n"); - for (int i = 0; i < 80; ++i) printf("="); printf("\n\n"); - + delete[] egvs; + puts("\n================================================================================\n"); return; } @@ -468,7 +489,7 @@ void DynMat::real2rec() for (int i = 0; i < 9; ++i) ibasevec[i] *= vol; - printf("\n"); for (int i = 0; i < 80; ++i) printf("="); + printf("\n================================================================================"); printf("\nBasis vectors of the unit cell in real space:"); for (int i = 0; i < sysdim; ++i){ printf("\n A%d: ", i+1); @@ -479,8 +500,7 @@ void DynMat::real2rec() printf("\n B%d: ", i+1); for (int j = 0; j < sysdim; ++j) printf("%8.4f ", ibasevec[i*3+j]); } - printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n"); - + puts("\n================================================================================"); return; } @@ -500,16 +520,17 @@ void DynMat::GaussJordan(int n, double *Mat) indxc = new int[n]; indxr = new int[n]; ipiv = new int[n]; - + + irow = icol = -1; for (i = 0; i < n; ++i) ipiv[i] = 0; for (i = 0; i < n; ++i){ - big = 0.; + big = 0.0; for (j = 0; j < n; ++j){ if (ipiv[j] != 1){ for (k = 0; k < n; ++k){ if (ipiv[k] == 0){ idr = j * n + k; - nmjk = abs(Mat[idr]); + nmjk = fabs(Mat[idr]); if (nmjk >= big){ big = nmjk; irow = j; @@ -602,6 +623,9 @@ void DynMat::help() printf(" will also inform the code to skip all q-points that is in the vicinity\n"); printf(" of the gamma point when evaluating phonon DOS and/or phonon dispersion.\n\n"); printf(" By default, this is not set; and not expected for uncharged systems.\n\n"); + printf(" -p prec To define the precision for symmetry identification with spglib.\n"); + printf(" By default, 1.e-3.\n\n"); + printf(" -save To record user input in `script.inp`, facilitating scripting.\n\n"); printf(" -h To print out this help info.\n\n"); printf(" file To define the filename that carries the binary dynamical matrice generated\n"); printf(" by fix-phonon. If not provided, the code will ask for it.\n"); @@ -680,13 +704,13 @@ void DynMat::Define_Conversion_Factor() * ---------------------------------------------------------------------------- */ void DynMat::ShowInfo() { - printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n"); + puts("\n================================================================================"); printf("Dynamical matrix is read from file: %s\n", binfile); printf("The system size in three dimension: %d x %d x %d\n", nx, ny, nz); printf("Number of atoms per unit cell : %d\n", nucell); printf("System dimension : %d\n", sysdim); printf("Boltzmann constant in used units : %g\n", boltz); - for (int i = 0; i < 80; ++i) printf("="); printf("\n"); + puts("================================================================================"); return; } /* --------------------------------------------------------------------*/ diff --git a/tools/phonon/dynmat.h b/tools/phonon/dynmat.h index 10caa1ce2a..e11badfde6 100644 --- a/tools/phonon/dynmat.h +++ b/tools/phonon/dynmat.h @@ -1,11 +1,9 @@ #ifndef DYNMAT_H #define DYNMAT_H -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "memory.h" -#include "interpolate.h" +#include "zheevd.h" + +#include class DynMat { public: @@ -15,7 +13,7 @@ public: int nx, ny, nz, nucell; int sysdim, fftdim; - double eml2f, eml2fc; + double eml2f, eml2fc, symprec; char *funit; void getDMq(double *); @@ -34,18 +32,18 @@ public: double **basis; int *attyp; + class UserInput *input; + private: int flag_skip, flag_reset_gamma; - Interpolate *interpolate; - - Memory *memory; + class Interpolate *interpolate; + class Memory *memory; - int nasr; void EnforceASR(); char *binfile, *dmfile; - double boltz, q[3]; + double boltz; doublecomplex **DM_all; diff --git a/tools/phonon/global.h b/tools/phonon/global.h index c20a91f317..bd1608ee07 100644 --- a/tools/phonon/global.h +++ b/tools/phonon/global.h @@ -1,7 +1,7 @@ #ifndef GLOBAL_H #define GLOBAL_H -#define ZERO 1.e-8 +#define ZERO 1.0e-8 #define MAXLINE 512 #define MIN(a,b) ((a)>(b)?(b):(a)) diff --git a/tools/phonon/green.cpp b/tools/phonon/green.cpp index 41f5b14886..0f687fd834 100644 --- a/tools/phonon/green.cpp +++ b/tools/phonon/green.cpp @@ -1,10 +1,10 @@ -#include -#include -#include -#include #include "green.h" + +#include "memory.h" + #include -#include "global.h" +#include +#include /******************************************************************************* * The class of Green is designed to evaluate the LDOS via the Green's Function @@ -59,7 +59,6 @@ Green::Green(const int ntm, const int sdim, const int niter, const double min, c dw = (wmax - wmin)/double(nw-1); memory->create(alpha, sysdim,nit, "Green_Green:alpha"); memory->create(beta, sysdim,nit+1,"Green_Green:beta"); - //memory->create(ldos, nw,sysdim, "Green_Green:ldos"); // use Lanczos algorithm to diagonalize the Hessian Lanczos(); @@ -224,8 +223,6 @@ void Green::recursion() { // local variables std::complex Z, rec_x, rec_x_inv; - std::complex cunit = std::complex(0.,1.); - double w = wmin; for (int i = 0; i < nw; ++i){ diff --git a/tools/phonon/green.h b/tools/phonon/green.h index 1a137e53ba..21e3b091dc 100644 --- a/tools/phonon/green.h +++ b/tools/phonon/green.h @@ -1,8 +1,6 @@ #ifndef GREEN_H #define GREEN_H -#include "memory.h" - class Green{ public: Green(const int, const int, const int, const double, const double, @@ -14,12 +12,11 @@ private: void Recursion(); void recursion(); - int ndos; double **ldos; int natom, iatom, sysdim, nit, nw, ndim; double dw, wmin, wmax, epson; double **alpha, **beta, **H; - Memory *memory; + class Memory *memory; }; #endif diff --git a/tools/phonon/input.cpp b/tools/phonon/input.cpp new file mode 100644 index 0000000000..c2059043c7 --- /dev/null +++ b/tools/phonon/input.cpp @@ -0,0 +1,35 @@ +#include "input.h" + +#include "global.h" + +/* ------------------------------------------------------------------- + * Constructor. If flag = 1, output user inputs as script.inp + * ---------------------------------------------------------------- */ +UserInput::UserInput(int flag) +{ + fp = NULL; + if (flag) fp = fopen("script.inp", "w"); + + return; +} + +/* ------------------------------------------------------------------- + * Deconstructor. Output user inputs as required and clear workspace. + * ---------------------------------------------------------------- */ +UserInput::~UserInput() +{ + if (fp) fclose(fp); + fp = NULL; +} + +/* ------------------------------------------------------------------- + * Read stdin and keep a record of it. + * ---------------------------------------------------------------- */ +void UserInput::read_stdin(char *str) +{ + fgets(str, MAXLINE, stdin); + if (fp) fprintf(fp, "%s", str); + + return; +} +/* ---------------------------------------------------------------- */ diff --git a/tools/phonon/input.h b/tools/phonon/input.h new file mode 100644 index 0000000000..c931895c4d --- /dev/null +++ b/tools/phonon/input.h @@ -0,0 +1,17 @@ +#ifndef INPUT_H +#define INPUT_H + +#include + +class UserInput { +public: + UserInput(int); + ~UserInput(); + + void read_stdin(char *); + +private: + FILE *fp; + +}; +#endif diff --git a/tools/phonon/interpolate.cpp b/tools/phonon/interpolate.cpp index 09e261c763..8ea551d1a8 100644 --- a/tools/phonon/interpolate.cpp +++ b/tools/phonon/interpolate.cpp @@ -1,6 +1,14 @@ + #include "interpolate.h" -#include "math.h" + #include "global.h" +#include "input.h" +#include "memory.h" +#include "tricubic.h" + +#include +#include +#include /* ---------------------------------------------------------------------------- * Constructor used to get info from caller, and prepare other necessary data @@ -19,6 +27,7 @@ Interpolate::Interpolate(int nx, int ny, int nz, int ndm, doublecomplex **DM) data = DM; Dfdx = Dfdy = Dfdz = D2fdxdy = D2fdxdz = D2fdydz = D3fdxdydz = NULL; flag_reset_gamma = flag_allocated_dfs = 0; + input = NULL; return; } @@ -265,17 +274,19 @@ void Interpolate::set_method() { char str[MAXLINE]; int im = 1; - printf("\n");for(int i=0; i<80; i++) printf("="); - printf("\nWhich interpolation method would you like to use?\n"); + if (input == NULL) input = new UserInput(0); + + puts("\n================================================================================"); + printf("Which interpolation method would you like to use?\n"); printf(" 1. Tricubic;\n 2. Trilinear;\n"); printf("Your choice [1]: "); - fgets(str,MAXLINE,stdin); + input->read_stdin(str); char *ptr = strtok(str," \t\n\r\f"); if (ptr) im = atoi(ptr); which =2-im%2; printf("Your selection: %d\n", which); - for(int i=0; i<80; i++) printf("="); printf("\n\n"); + puts("================================================================================\n"); if (which == 1) tricubic_init(); diff --git a/tools/phonon/interpolate.h b/tools/phonon/interpolate.h index c650b30908..b9e0242b96 100644 --- a/tools/phonon/interpolate.h +++ b/tools/phonon/interpolate.h @@ -1,16 +1,7 @@ #ifndef INTERPOLATION_H #define INTERPOLATION_H -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "memory.h" -#include "tricubic.h" - -extern "C"{ -#include "f2c.h" -#include "clapack.h" -} +#include "zheevd.h" class Interpolate{ public: @@ -23,11 +14,13 @@ public: int UseGamma; + class UserInput *input; + private: void tricubic_init(); void tricubic(double *, doublecomplex *); void trilinear(double *, doublecomplex *); - Memory *memory; + class Memory *memory; int which; int Nx, Ny, Nz, Npt, ndim; diff --git a/tools/phonon/kpath.cpp b/tools/phonon/kpath.cpp index 842e680611..49730b42b6 100644 --- a/tools/phonon/kpath.cpp +++ b/tools/phonon/kpath.cpp @@ -1,11 +1,19 @@ -#include "global.h" + #include "kpath.h" +#include "global.h" +#include "dynmat.h" +#include "memory.h" +#include "qnodes.h" + #ifdef UseSPG extern "C"{ #include "spglib.h" } -#include "math.h" +#include +#include +#include +#include /* ---------------------------------------------------------------------------- * Class kPath will help to find the high symmetry k-path for a given lattice. @@ -47,11 +55,15 @@ kPath::kPath(DynMat *dm, QNodes *qn) for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim]; // get the space group number - double symprec = 1.e-4, pos[num_atom][3]; + double symprec = 1.0e-3; + double **pos; + memory->create(pos,num_atom,3,"kpath:pos"); + if (dynmat->symprec > 0.0) symprec = dynmat->symprec; + for (int i = 0; i < num_atom; ++i) for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j]; - spgnum = spg_get_international(symbol, latvec, pos, attyp, num_atom, symprec); - + spgnum = spg_get_international(symbol, latvec, (double (*)[3])pos, attyp, num_atom, symprec); + memory->destroy(pos); return; } @@ -61,7 +73,7 @@ kPath::kPath(DynMat *dm, QNodes *qn) void kPath::show_info() { // display the unit cell info read - for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n"); + puts("--------------------------------------------------------------------------------"); printf("The basis vectors of the unit cell:\n"); for (int idim = 0; idim < 3; ++idim){ printf(" A%d =", idim+1); @@ -76,12 +88,10 @@ void kPath::show_info() if (num_atom > NUMATOM) printf(" ... (%d atoms omitted.)\n", num_atom-NUMATOM); printf("The space group number of your unit cell is: %d => %s\n", spgnum, symbol); - for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n"); - + puts("--------------------------------------------------------------------------------"); return; } - /* ---------------------------------------------------------------------------- * Free the memeory used by kPath. * ---------------------------------------------------------------------------- */ @@ -2765,9 +2775,16 @@ void kPath::show_path() if (q == NULL) return; int nbin = q->ndstr.size(); if (nbin > 0){ - printf("\nk-path for the current lattice will be:\n\t%s", q->ndstr[0].c_str()); + puts("\n--------------------------------------------------------------------------------"); + printf("k-path for the current lattice will be:\n %s", q->ndstr[0].c_str()); for (int is = 1; is < nbin; ++is) printf("-%s", q->ndstr[is].c_str()); - printf("\n"); + + printf("\n\nThe fractional coordinates of these paths are:\n"); + for (int is = 0; is < nbin-1; ++is) + printf(" [%6.4f %6.4f %6.4f] --> [%6.4f %6.4f %6.4f] (%s - %s)\n", q->qs[is][0], + q->qs[is][1], q->qs[is][2], q->qe[is][0], q->qe[is][1], q->qe[is][2], + q->ndstr[is].c_str(), q->ndstr[is+1].c_str() ); + puts("--------------------------------------------------------------------------------"); } return; diff --git a/tools/phonon/kpath.h b/tools/phonon/kpath.h index bbcedf15d4..ade1e3a27c 100644 --- a/tools/phonon/kpath.h +++ b/tools/phonon/kpath.h @@ -4,14 +4,9 @@ #ifndef KPATH_H #define KPATH_H -#include "qnodes.h" -#include "dynmat.h" -#include "memory.h" - class kPath{ public: - - kPath(DynMat *, QNodes *); + kPath(class DynMat *, class QNodes *); ~kPath(); void kpath(); @@ -19,13 +14,11 @@ public: void show_info(); private: - - Memory *memory; - - DynMat *dynmat; - QNodes *q; + class Memory *memory; + class DynMat *dynmat; + class QNodes *q; char symbol[11]; - int spgnum, sysdim, fftdim, num_atom, *attyp; + int spgnum, sysdim, num_atom, *attyp; double latvec[3][3], **atpos; }; diff --git a/tools/phonon/main.cpp b/tools/phonon/main.cpp index d7d5baa1cc..d0193a037c 100644 --- a/tools/phonon/main.cpp +++ b/tools/phonon/main.cpp @@ -1,10 +1,6 @@ -#include "stdio.h" -#include "stdlib.h" #include "dynmat.h" #include "phonon.h" -using namespace std; - int main(int argc, char** argv) { diff --git a/tools/phonon/memory.cpp b/tools/phonon/memory.cpp index 4d65e83ef8..18dfa2fa43 100644 --- a/tools/phonon/memory.cpp +++ b/tools/phonon/memory.cpp @@ -1,8 +1,8 @@ -#include "stdio.h" -#include "stdlib.h" -#include "string.h" #include "memory.h" +#include +#include + /* ---------------------------------------------------------------------- safe malloc ------------------------------------------------------------------------- */ diff --git a/tools/phonon/memory.h b/tools/phonon/memory.h index ae2feceba3..13eeca4b14 100644 --- a/tools/phonon/memory.h +++ b/tools/phonon/memory.h @@ -4,18 +4,16 @@ #define __STDC_LIMIT_MACROS #define __STDC_FORMAT_MACROS -#include "stdio.h" -#include "stdlib.h" -#include "limits.h" -#include "stdint.h" -#include "inttypes.h" +#include +#include +#include typedef int64_t bigint; #define BIGINT_FORMAT "%" PRId64 #define ATOBIGINT atoll class Memory { - public: +public: Memory(){}; void *smalloc(bigint n, const char *); @@ -24,11 +22,11 @@ class Memory { void fail(const char *); /* ---------------------------------------------------------------------- - create a 1d array -------------------------------------------------------------------------- */ + create a 1d array + ------------------------------------------------------------------------- */ template - TYPE *create(TYPE *&array, int n, const char *name) + TYPE *create(TYPE *&array, int n, const char *name) { bigint nbytes = sizeof(TYPE) * n; array = (TYPE *) smalloc(nbytes,name); @@ -36,42 +34,42 @@ class Memory { }; template - TYPE **create(TYPE **&array, int n, const char *name) {fail(name);} + TYPE **create(TYPE **&, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- grow or shrink 1d array -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - TYPE *grow(TYPE *&array, int n, const char *name) + TYPE *grow(TYPE *&array, int n, const char *name) { if (array == NULL) return create(array,n,name); - + bigint nbytes = sizeof(TYPE) * n; array = (TYPE *) srealloc(array,nbytes,name); return array; }; template - TYPE **grow(TYPE **&array, int n, const char *name) {fail(name);} + TYPE **grow(TYPE **&, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- - destroy a 1d array -------------------------------------------------------------------------- */ + destroy a 1d array + ------------------------------------------------------------------------- */ template - void destroy(TYPE *array) + void destroy(TYPE *array) { sfree(array); }; /* ---------------------------------------------------------------------- - create a 1d array with index from nlo to nhi inclusive + create a 1d array with index from nlo to nhi inclusive cannot grow it -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name) + TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name) { bigint nbytes = sizeof(TYPE) * (nhi-nlo+1); array = (TYPE *) smalloc(nbytes,name); @@ -80,76 +78,76 @@ class Memory { } template - TYPE **create1d_offset(TYPE **&array, int nlo, int nhi, const char *name) + TYPE **create1d_offset(TYPE **&, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- - destroy a 1d array with index offset -------------------------------------------------------------------------- */ + destroy a 1d array with index offset + ------------------------------------------------------------------------- */ template - void destroy1d_offset(TYPE *array, int offset) + void destroy1d_offset(TYPE *array, int offset) { if (array) sfree(&array[offset]); } /* ---------------------------------------------------------------------- - create a 2d array -------------------------------------------------------------------------- */ + create a 2d array + ------------------------------------------------------------------------- */ template - TYPE **create(TYPE **&array, int n1, int n2, const char *name) + TYPE **create(TYPE **&array, int n1, int n2, const char *name) { bigint nbytes = sizeof(TYPE) * n1*n2; TYPE *data = (TYPE *) smalloc(nbytes,name); nbytes = sizeof(TYPE *) * n1; array = (TYPE **) smalloc(nbytes,name); - + int n = 0; for (int i = 0; i < n1; i++) { - array[i] = &data[n]; - n += n2; + array[i] = &data[n]; + n += n2; } return array; } template - TYPE ***create(TYPE ***&array, int n1, int n2, const char *name) + TYPE ***create(TYPE ***&, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- grow or shrink 1st dim of a 2d array last dim must stay the same -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - TYPE **grow(TYPE **&array, int n1, int n2, const char *name) + TYPE **grow(TYPE **&array, int n1, int n2, const char *name) { if (array == NULL) return create(array,n1,n2,name); - + bigint nbytes = sizeof(TYPE) * n1*n2; TYPE *data = (TYPE *) srealloc(array[0],nbytes,name); nbytes = sizeof(TYPE *) * n1; array = (TYPE **) srealloc(array,nbytes,name); - + int n = 0; for (int i = 0; i < n1; i++) { - array[i] = &data[n]; - n += n2; + array[i] = &data[n]; + n += n2; } return array; } template - TYPE ***grow(TYPE ***&array, int n1, int n2, const char *name) + TYPE ***grow(TYPE ***&, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- - destroy a 2d array -------------------------------------------------------------------------- */ + destroy a 2d array + ------------------------------------------------------------------------- */ template - void destroy(TYPE **array) + void destroy(TYPE **array) { if (array == NULL) return; sfree(array[0]); @@ -157,42 +155,11 @@ class Memory { } /* ---------------------------------------------------------------------- - create a 2d array with 2nd index from n2lo to n2hi inclusive - cannot grow it -------------------------------------------------------------------------- */ + create a 3d array + ------------------------------------------------------------------------- */ template - TYPE **create2d_offset(TYPE **&array, int n1, int n2lo, int n2hi, - const char *name) - { - int n2 = n2hi - n2lo + 1; - create(array,n1,n2,name); - for (int i = 0; i < n1; i++) array[i] -= n2lo; - return array; - } - - template - TYPE ***create2d_offset(TYPE ***&array, int n1, int n2lo, int n2hi, - const char *name) {fail(name);} - -/* ---------------------------------------------------------------------- - destroy a 2d array with 2nd index offset -------------------------------------------------------------------------- */ - - template - void destroy2d_offset(TYPE **array, int offset) - { - if (array == NULL) return; - sfree(&array[0][offset]); - sfree(array); - } - -/* ---------------------------------------------------------------------- - create a 3d array -------------------------------------------------------------------------- */ - - template - TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name) + TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name) { bigint nbytes = sizeof(TYPE) * n1*n2*n3; TYPE *data = (TYPE *) smalloc(nbytes,name); @@ -200,62 +167,62 @@ class Memory { TYPE **plane = (TYPE **) smalloc(nbytes,name); nbytes = sizeof(TYPE **) * n1; array = (TYPE ***) smalloc(nbytes,name); - + int i,j; int n = 0; for (i = 0; i < n1; i++) { - array[i] = &plane[i*n2]; - for (j = 0; j < n2; j++) { - plane[i*n2+j] = &data[n]; - n += n3; - } + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &data[n]; + n += n3; + } } return array; } template - TYPE ****create(TYPE ****&array, int n1, int n2, int n3, const char *name) + TYPE ****create(TYPE ****&, int, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- grow or shrink 1st dim of a 3d array last 2 dims must stay the same -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name) + TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name) { if (array == NULL) return create(array,n1,n2,n3,name); - + bigint nbytes = sizeof(TYPE) * n1*n2*n3; TYPE *data = (TYPE *) srealloc(array[0][0],nbytes,name); nbytes = sizeof(TYPE *) * n1*n2; TYPE **plane = (TYPE **) srealloc(array[0],nbytes,name); nbytes = sizeof(TYPE **) * n1; array = (TYPE ***) srealloc(array,nbytes,name); - + int i,j; int n = 0; for (i = 0; i < n1; i++) { - array[i] = &plane[i*n2]; - for (j = 0; j < n2; j++) { - plane[i*n2+j] = &data[n]; - n += n3; - } + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &data[n]; + n += n3; + } } return array; } template - TYPE ****grow(TYPE ****&array, int n1, int n2, int n3, const char *name) + TYPE ****grow(TYPE ****&, int, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- - destroy a 3d array -------------------------------------------------------------------------- */ + destroy a 3d array + ------------------------------------------------------------------------- */ template - void destroy(TYPE ***array) + void destroy(TYPE ***array) { if (array == NULL) return; sfree(array[0][0]); @@ -263,168 +230,23 @@ class Memory { sfree(array); } -/* ---------------------------------------------------------------------- - create a 3d array with 1st index from n1lo to n1hi inclusive - cannot grow it -------------------------------------------------------------------------- */ - - template - TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, - int n2, int n3, const char *name) - { - int n1 = n1hi - n1lo + 1; - create(array,n1,n2,n3,name); - array -= n1lo; - return array; - } - - template - TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi, - int n2, int n3, const char *name) - {fail(name);} - -/* ---------------------------------------------------------------------- - free a 3d array with 1st index offset -------------------------------------------------------------------------- */ - - template - void destroy3d_offset(TYPE ***array, int offset) - { - if (array) destroy(&array[offset]); - } - -/* ---------------------------------------------------------------------- - create a 3d array with - 1st index from n1lo to n1hi inclusive, - 2nd index from n2lo to n2hi inclusive, - 3rd index from n3lo to n3hi inclusive - cannot grow it -------------------------------------------------------------------------- */ - - template - TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, - int n2lo, int n2hi, int n3lo, int n3hi, - const char *name) - { - int n1 = n1hi - n1lo + 1; - int n2 = n2hi - n2lo + 1; - int n3 = n3hi - n3lo + 1; - create(array,n1,n2,n3,name); - - for (int i = 0; i < n1*n2; i++) array[0][i] -= n3lo; - for (int i = 0; i < n1; i++) array[i] -= n2lo; - array -= n1lo; - return array; - } - - template - TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi, - int n2lo, int n2hi, int n3lo, int n3hi, - const char *name) - {fail(name);} - -/* ---------------------------------------------------------------------- - free a 3d array with all 3 indices offset -------------------------------------------------------------------------- */ - - template - void destroy3d_offset(TYPE ***array, - int n1_offset, int n2_offset, int n3_offset) - { - if (array == NULL) return; - sfree(&array[n1_offset][n2_offset][n3_offset]); - sfree(&array[n1_offset][n2_offset]); - sfree(&array[n1_offset]); - } - -/* ---------------------------------------------------------------------- - create a 4d array -------------------------------------------------------------------------- */ - - template - TYPE ****create(TYPE ****&array, int n1, int n2, int n3, int n4, - const char *name) - { - bigint nbytes = sizeof(TYPE) * n1*n2*n3*n4; - TYPE *data = (double *) smalloc(nbytes,name); - nbytes = sizeof(TYPE *) * n1*n2*n3; - TYPE **cube = (double **) smalloc(nbytes,name); - nbytes = sizeof(TYPE **) * n1*n2; - TYPE ***plane = (double ***) smalloc(nbytes,name); - nbytes = sizeof(TYPE ***) * n1; - array = (double ****) smalloc(nbytes,name); - - int i,j,k; - int n = 0; - for (i = 0; i < n1; i++) { - array[i] = &plane[i*n2]; - for (j = 0; j < n2; j++) { - plane[i*n2+j] = &cube[i*n2*n3+j*n3]; - for (k = 0; k < n3; k++) { - cube[i*n2*n3+j*n3+k] = &data[n]; - n += n4; - } - } - } - return array; - } - - template - TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4, - const char *name) - {fail(name);} - -/* ---------------------------------------------------------------------- - destroy a 4d array -------------------------------------------------------------------------- */ - - template - void destroy(TYPE ****array) - { - if (array == NULL) return; - sfree(array[0][0][0]); - sfree(array[0][0]); - sfree(array[0]); - sfree(array); - } - /* ---------------------------------------------------------------------- memory usage of arrays, including pointers -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - bigint usage(TYPE *array, int n) + bigint usage(TYPE *, int n) { bigint bytes = sizeof(TYPE) * n; return bytes; } template - bigint usage(TYPE **array, int n1, int n2) + bigint usage(TYPE **, int n1, int n2) { bigint bytes = sizeof(TYPE) * n1*n2; bytes += sizeof(TYPE *) * n1; return bytes; } - - template - bigint usage(TYPE ***array, int n1, int n2, int n3) - { - bigint bytes = sizeof(TYPE) * n1*n2*n3; - bytes += sizeof(TYPE *) * n1*n2; - bytes += sizeof(TYPE **) * n1; - return bytes; - } - - template - bigint usage(TYPE ****array, int n1, int n2, int n3, int n4) - { - bigint bytes = sizeof(TYPE) * n1*n2*n3*n4; - bytes += sizeof(TYPE *) * n1*n2*n3; - bytes += sizeof(TYPE **) * n1*n2; - bytes += sizeof(TYPE ***) * n1; - return bytes; - } }; - #endif diff --git a/tools/phonon/phonon.cpp b/tools/phonon/phonon.cpp index 56fa409e06..06372dcd1b 100644 --- a/tools/phonon/phonon.cpp +++ b/tools/phonon/phonon.cpp @@ -1,9 +1,18 @@ -#include -#include "string.h" + #include "phonon.h" -#include "green.h" -#include "timer.h" + #include "global.h" +#include "dynmat.h" +#include "green.h" +#include "input.h" +#include "memory.h" +#include "timer.h" +#include "zheevd.h" + +#include +#include +#include +#include #ifdef UseSPG extern "C"{ @@ -27,6 +36,7 @@ Phonon::Phonon(DynMat *dm) dynmat = dm; sysdim = dynmat->sysdim; ndim = dynmat->fftdim; + input = dm->input; dos = NULL; ldos = NULL; qpts = NULL; @@ -42,10 +52,7 @@ Phonon::Phonon(DynMat *dm) // display the menu char str[MAXLINE]; while ( 1 ){ - printf("\n"); - for (int i = 0; i < 37; ++i) printf("="); - printf(" Menu "); - for (int i = 0; i < 37; ++i) printf("="); printf("\n"); + puts("\n===================================== Menu ====================================="); printf(" 1. Phonon DOS evaluation;\n"); printf(" 2. Phonon dispersion curves;\n"); printf(" 3. Dynamical matrix at arbitrary q;\n"); @@ -64,9 +71,10 @@ Phonon::Phonon(DynMat *dm) // read user choice int job = 0; printf("Your choice [0]: "); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) job = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) job = atoi(strtok(str," \t\n\r\f")); printf("\nYour selection: %d\n", job); - for (int i = 0; i < 80; ++i) printf("=");printf("\n\n"); + puts("================================================================================\n"); // now to do the job according to user's choice if (job == 1) pdos(); @@ -138,7 +146,8 @@ void Phonon::pdos() // Now to ask for the output frequency range printf("\nThe frequency range of all q-points are: [%g %g]\n", fmin, fmax); printf("Please input the desired range to get DOS [%g %g]: ", fmin, fmax); - if (count_words(fgets(str,MAXLINE,stdin)) >= 2){ + input->read_stdin(str); + if (count_words(str) >= 2){ fmin = atof(strtok(str," \t\n\r\f")); fmax = atof(strtok(NULL," \t\n\r\f")); } @@ -147,7 +156,8 @@ void Phonon::pdos() ndos = 201; printf("Please input the number of intervals [%d]: ", ndos); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f")); ndos += (ndos+1)%2; ndos = MAX(2,ndos); @@ -170,7 +180,8 @@ void Phonon::pdos() // smooth dos ? printf("Would you like to smooth the phonon dos? (y/n)[n]: "); - if (count_words(fgets(str,MAXLINE,stdin)) > 0){ + input->read_stdin(str); + if (count_words(str) > 0){ char *flag = strtok(str," \t\n\r\f"); if (strcmp(flag,"y") == 0 || strcmp(flag,"Y") == 0) smooth(dos, ndos); } @@ -194,7 +205,8 @@ void Phonon::writeDOS() char str[MAXLINE]; // now to output the phonon DOS printf("\nPlease input the filename to write DOS [pdos.dat]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdos.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "pdos.dat"); char *fname = strtok(str," \t\n\r\f"); printf("The total phonon DOS will be written to file: %s\n", fname); @@ -234,7 +246,7 @@ void Phonon::writeLDOS() const double one3 = 1./double(sysdim); char str[MAXLINE]; for (int ilocal = 0; ilocal < nlocal; ++ilocal){ - sprintf(str,"pldos_%d.dat", locals[ilocal]); + snprintf(str, MAXLINE-1, "pldos_%d.dat", locals[ilocal]); char *fname = strtok(str," \t\n\r\f"); FILE *fp = fopen(fname, "w"); fname = NULL; @@ -281,7 +293,7 @@ void Phonon::ldos_rsgf() fmin = fmax = egvs[0]; for (int i = 1; i < ndim; ++i){fmin = MIN(fmin, egvs[i]); fmax = MAX(fmax, egvs[i]);} - delete []egvs; + delete[] egvs; } else { fmin = 0.; fmax = 20.; @@ -297,7 +309,8 @@ void Phonon::ldos_rsgf() printf("\nThere are %d atoms in each unit cell of your lattice.\n", dynmat->nucell); printf("Please input the index/index range/index range and increment of atom(s)\n"); printf("in the unit cell to evaluate LDOS, q to exit [%d]: ", ik); - int nr = count_words( fgets(str,MAXLINE,stdin) ); + input->read_stdin(str); + int nr = count_words(str); if (nr < 1){ istr = iend = ik; iinc = 1; @@ -327,7 +340,8 @@ void Phonon::ldos_rsgf() } printf("Please input the frequency range to evaluate LDOS [%g %g]: ", fmin, fmax); - if (count_words(fgets(str,MAXLINE,stdin)) >= 2){ + input->read_stdin(str); + if (count_words(str) >= 2){ fmin = atof(strtok(str," \t\n\r\f")); fmax = atof(strtok(NULL," \t\n\r\f")); } @@ -335,16 +349,19 @@ void Phonon::ldos_rsgf() printf("The frequency range for your LDOS is [%g %g].\n", fmin, fmax); printf("Please input the desired number of points in LDOS [%d]: ", ndos); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f")); if (ndos < 2) break; ndos += (ndos+1)%2; printf("Please input the maximum # of Lanczos iterations [%d]: ", nit); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) nit = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) nit = atoi(strtok(str," \t\n\r\f")); if (nit < 1) break; printf("Please input the value of epsilon for delta-function [%g]: ", eps); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) eps = atof(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) eps = atof(strtok(str," \t\n\r\f")); if (eps <= 0.) break; // prepare array for local pdos @@ -395,8 +412,11 @@ void Phonon::dmanyq() { char str[MAXLINE]; double q[3]; - do printf("Please input the q-point to output the dynamical matrix:"); - while (count_words(fgets(str,MAXLINE,stdin)) < 3); + while ( 1 ){ + printf("Please input the q-point to output the dynamical matrix: "); + input->read_stdin(str); + if (count_words(str) >= 3) break; + } q[0] = atof(strtok(str," \t\n\r\f")); q[1] = atof(strtok(NULL," \t\n\r\f")); q[2] = atof(strtok(NULL," \t\n\r\f")); @@ -413,11 +433,13 @@ void Phonon::dmanyq() void Phonon::vfanyq() { char str[MAXLINE]; - double q[3], egvs[ndim]; + double q[3]; + double *egvs = new double[ndim]; while ( 1 ){ printf("Please input the q-point to compute the frequencies, q to exit: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 3) break; + input->read_stdin(str); + if (count_words(str) < 3) break; q[0] = atof(strtok(str, " \t\n\r\f")); q[1] = atof(strtok(NULL," \t\n\r\f")); @@ -427,9 +449,11 @@ void Phonon::vfanyq() dynmat->geteigen(egvs, 0); printf("q-point: [%lg %lg %lg], ", q[0], q[1], q[2]); printf("vibrational frequencies at this q-point:\n"); - for (int i = 0; i < ndim; ++i) printf("%lg ", egvs[i]); printf("\n\n"); + for (int i = 0; i < ndim; ++i) printf("%lg ", egvs[i]); + printf("\n\n"); } - + + delete[] egvs; return; } @@ -439,15 +463,18 @@ void Phonon::vfanyq() void Phonon::vecanyq() { char str[MAXLINE]; - double q[3], egvs[ndim]; + double q[3]; + double *egvs = new double[ndim]; doublecomplex **eigvec = dynmat->DM_q; printf("Please input the filename to output the result [eigvec.dat]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str,"eigvec.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str,"eigvec.dat"); FILE *fp = fopen(strtok(str," \t\n\r\f"), "w"); while ( 1 ){ printf("Please input the q-point to compute the frequencies, q to exit: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 3) break; + input->read_stdin(str); + if (count_words(str) < 3) break; q[0] = atof(strtok(str, " \t\n\r\f")); q[1] = atof(strtok(NULL," \t\n\r\f")); @@ -475,6 +502,7 @@ void Phonon::vecanyq() fprintf(fp,"\n"); } fclose(fp); + delete[] egvs; eigvec = NULL; return; } @@ -488,7 +516,8 @@ void Phonon::DMdisp() char str[MAXLINE]; printf("Please input the filename to output the DM data [DMDisp.dat]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "DMDisp.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "DMDisp.dat"); char *fname = strtok(str," \t\n\r\f"); FILE *fp = fopen(fname, "w"); fname = NULL; @@ -503,7 +532,8 @@ void Phonon::DMdisp() for (int i = 0; i < 3; ++i) qstr[i] = qend[i]; printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]); - int n = count_words(fgets(str,MAXLINE,stdin)); + input->read_stdin(str); + int n = count_words(str); char *ptr = strtok(str," \t\n\r\f"); if ((n == 1) && (strcmp(ptr,"q") == 0)) break; else if (n >= 3){ @@ -512,14 +542,18 @@ void Phonon::DMdisp() qstr[2] = atof(strtok(NULL," \t\n\r\f")); } - do printf("Please input the end q-point in unit of B1->B3: "); - while (count_words(fgets(str,MAXLINE,stdin)) < 3); + while ( 1 ){ + printf("Please input the end q-point in unit of B1->B3: "); + input->read_stdin(str); + if (count_words(str) >= 3) break; + } qend[0] = atof(strtok(str," \t\n\r\f")); qend[1] = atof(strtok(NULL," \t\n\r\f")); qend[2] = atof(strtok(NULL," \t\n\r\f")); printf("Please input the # of points along the line [%d]: ", nq); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) nq = atoi(strtok(str," \t\n\r\f")); nq = MAX(nq,2); for (int i=0; i<3; i++) qinc[i] = (qend[i]-qstr[i])/double(nq-1); @@ -588,7 +622,8 @@ void Phonon::therm() char str[MAXLINE]; printf("\nPlease input the filename to output thermal properties [therm.dat]:"); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "therm.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "therm.dat"); char *fname = strtok(str," \t\n\r\f"); FILE *fp = fopen(fname, "a"); fname = NULL; // header line @@ -630,7 +665,8 @@ void Phonon::therm() fprintf(fp,"%lg %lg %lg %lg %lg %lg\n", T, Uvib, Svib, Fvib, ZPE, Cvib); printf("Please input the desired temperature (K), enter to exit: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) break; + input->read_stdin(str); + if (count_words(str) < 1) break; T = atof(strtok(str," \t\n\r\f")); } while (T > 0.); @@ -646,12 +682,14 @@ void Phonon::local_therm() { char str[MAXLINE]; printf("\nWould you like to compute the local thermal properties (y/n)[n]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) return; + input->read_stdin(str); + if (count_words(str) < 1) return; char *ptr = strtok(str," \t\n\r\f"); if (strcmp(ptr,"y") != 0 && strcmp(ptr, "Y") != 0 && strcmp(ptr, "yes") != 0) return; printf("Please input the filename to output vibrational thermal info [localtherm.dat]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "localtherm.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "localtherm.dat"); FILE *fp = fopen(strtok(str," \t\n\r\f"), "w"); fprintf(fp,"# atom Temp U_vib (eV) S_vib (kB) F_vib (eV) C_vib (kB) ZPE (eV)\n"); @@ -672,7 +710,8 @@ void Phonon::local_therm() while ( 1 ){ printf("\nPlease input the temperature at which to evaluate the local vibrational\n"); printf("thermal properties, non-positive number to exit [%g]: ", T); - if (count_words(fgets(str,MAXLINE,stdin)) > 0){ + input->read_stdin(str); + if (count_words(str) > 0){ T = atoi(strtok(str," \t\n\r\f")); if (T <= 0.) break; } @@ -765,7 +804,8 @@ void Phonon::QMesh() printf("\nThe q-mesh size from the read dynamical matrix is: %d x %d x %d\n", nx, ny, nz); printf("A denser mesh can be interpolated, but NOTE a too dense mesh can cause segmentation fault.\n"); printf("Please input your desired q-mesh size [%d %d %d]: ", nx, ny, nz); - if (count_words(fgets(str,MAXLINE,stdin)) >= 3){ + input->read_stdin(str); + if (count_words(str) >= 3){ nx = atoi(strtok(str," \t\n\r\f")); ny = atoi(strtok(NULL," \t\n\r\f")); nz = atoi(strtok(NULL," \t\n\r\f")); @@ -780,7 +820,8 @@ void Phonon::QMesh() int method = 2; printf("Please select your method to generate the q-points:\n"); printf(" 1. uniform;\n 2. Monkhost-Pack mesh;\nYour choice [2]: "); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) method = atoi(strtok(str," \t\n\r\f")); method = 2 - method%2; printf("Your selection: %d\n", method); #endif @@ -831,7 +872,7 @@ void Phonon::QMesh() for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim]; // display the unit cell info read - printf("\n");for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + puts("\n================================================================================"); printf("The basis vectors of the unit cell:\n"); for (int idim = 0; idim < 3; ++idim) printf(" A%d = %lg %lg %lg\n", idim+1, latvec[0][idim], latvec[1][idim], latvec[2][idim]); @@ -845,14 +886,21 @@ void Phonon::QMesh() mesh[0] = nx; mesh[1] = ny; mesh[2] = nz; shift[0] = shift[1] = shift[2] = 0; int num_grid = mesh[0]*mesh[1]*mesh[2]; - int grid_point[num_grid][3], map[num_grid]; - double symprec = 1.e-4, pos[num_atom][3]; + int **grid_point; + memory->create(grid_point, num_grid, 3, "phonon:grid_point"); + int *map = new int[num_grid]; + double symprec = 1.0e-3; + double **pos; + memory->create(pos, num_atom, 3, "phonon:pos"); + if (dynmat->symprec > 0.) symprec = dynmat->symprec; for (int i = 0; i < num_atom; ++i) - for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j]; + for (int j = 0; j < 3; ++j) + pos[i][j] = atpos[i][j]; // if spglib >= 1.0.3 is used - nq = spg_get_ir_reciprocal_mesh(grid_point, map, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec); + nq = spg_get_ir_reciprocal_mesh((int (*)[3])grid_point, map, mesh, shift, is_time_reversal, + latvec, (double (*)[3])pos, attyp, num_atom, symprec); memory->create(wt, nq, "QMesh:wt"); memory->create(qpts, nq,3,"QMesh:qpts"); @@ -873,11 +921,14 @@ void Phonon::QMesh() qpts[numq][2] = double(grid_point[i][2])/double(mesh[2]); numq++; } - wt[iq2idx[iq]] += 1.; + wt[iq2idx[iq]] += 1.0; } - delete []iq2idx; - - double wsum = 0.; + delete[] iq2idx; + delete[] map; + memory->destroy(grid_point); + memory->destroy(pos); + + double wsum = 0.0; for (int iq = 0; iq < nq; ++iq) wsum += wt[iq]; for (int iq = 0; iq < nq; ++iq) wt[iq] /= wsum; @@ -898,7 +949,8 @@ void Phonon::ldos_egv() char str[MAXLINE], *ptr; printf("\nThe # of atoms per cell is: %d, please input the atom IDs to compute\n", dynmat->nucell); printf("local PDOS, IDs begin with 0: "); - int nmax = count_words(fgets(str,MAXLINE,stdin)); + input->read_stdin(str); + int nmax = count_words(str); if (nmax < 1) return; memory->destroy(locals); @@ -920,7 +972,8 @@ void Phonon::ldos_egv() fmin = 0.; fmax = 10.; printf("Please input the freqency (nv, THz) range to compute PDOS [%g %g]: ", fmin, fmax); - if (count_words(fgets(str,MAXLINE,stdin)) >= 2) { + input->read_stdin(str); + if (count_words(str) >= 2) { fmin = atof(strtok(str," \t\n\r\f")); fmax = atof(strtok(NULL," \t\n\r\f")); } @@ -928,7 +981,8 @@ void Phonon::ldos_egv() ndos = 201; printf("Please input your desired # of points in PDOS [%d]: ", ndos); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f")); if (ndos < 2) return; ndos += (ndos+1)%2; @@ -957,7 +1011,8 @@ void Phonon::ldos_egv() Timer *time = new Timer(); // memory and pointer for eigenvalues and eigenvectors - double egval[ndim], offset=fmin-0.5*df; + double offset=fmin-0.5*df; + double *egval = new double[ndim]; doublecomplex **egvec = dynmat->DM_q; printf("\nNow to compute the phonons and DOSs "); fflush(stdout); @@ -985,6 +1040,7 @@ void Phonon::ldos_egv() } } } + delete[] egval; egvec = NULL; printf("Done!\nNow to normalize the DOSs ..."); fflush(stdout); @@ -1008,10 +1064,7 @@ void Phonon::ldos_egv() * ---------------------------------------------------------------------------- */ void Phonon::ShowCell() { - printf("\n"); - for (int i = 0; i < 30; ++i) printf("="); - printf(" Unit Cell Info "); - for (int i = 0; i < 30; ++i) printf("="); printf("\n"); + puts("============================== Unit Cell Info =============================="); printf("Number of atoms in the unit cell: %d\n", dynmat->nucell); printf("Basis vectors of the unit cell:\n"); printf(" %15.8f %15.8f %15.8f\n", dynmat->basevec[0], dynmat->basevec[1], dynmat->basevec[2]); @@ -1024,8 +1077,7 @@ void Phonon::ShowCell() printf("Atomic type and fractional coordinates:\n"); for (int i = 0; i < dynmat->nucell; ++i) printf("%4d %12.8f %12.8f %12.8f\n", dynmat->attyp[i], dynmat->basis[i][0], dynmat->basis[i][1], dynmat->basis[i][2]); - for (int i = 0; i < 80; ++i) printf("="); - printf("\n"); + puts("================================================================================"); return; } @@ -1101,7 +1153,7 @@ int Phonon::count_words(const char *line) strcpy(copy,line); char *ptr; - if (ptr = strchr(copy,'#')) *ptr = '\0'; + if ((ptr = strchr(copy,'#'))) *ptr = '\0'; if (strtok(copy," \t\n\r\f") == NULL) { memory->destroy(copy); diff --git a/tools/phonon/phonon.h b/tools/phonon/phonon.h index 69a1fe5d50..46ce17f268 100644 --- a/tools/phonon/phonon.h +++ b/tools/phonon/phonon.h @@ -1,22 +1,16 @@ #ifndef PHONON_H #define PHONON_H -#include "stdio.h" -#include "stdlib.h" -#include -#include "dynmat.h" -#include "memory.h" - -using namespace std; - class Phonon{ public: - Phonon(DynMat *); + Phonon(class DynMat *); ~Phonon(); - DynMat *dynmat; + class DynMat *dynmat; private: + class UserInput *input; + int nq, ndim, sysdim; double **qpts, *wt; double **eigs; @@ -25,7 +19,7 @@ private: double *dos, fmin, fmax, df, rdf; double ***ldos; - Memory *memory; + class Memory *memory; void QMesh(); void ComputeAll(); diff --git a/tools/phonon/phonopy.cpp b/tools/phonon/phonopy.cpp index cfe32ab61a..a6b991efac 100644 --- a/tools/phonon/phonopy.cpp +++ b/tools/phonon/phonopy.cpp @@ -1,9 +1,20 @@ + #ifdef FFTW3 -#include + #include "phonopy.h" -#include "math.h" -#include "kpath.h" -#include "fftw3.h" + +#include "global.h" +#include "dynmat.h" +#include "input.h" +#include "memory.h" + +#include + +#include +#include +#include +#include +#include /* ---------------------------------------------------------------------------- * Class Phonopy is designed to interface with phonopy. @@ -14,13 +25,16 @@ Phonopy::Phonopy(DynMat *dynmat) memory = new Memory(); sysdim = dm->sysdim; fftdim = dm->fftdim; + input = dm->input; fftdim2 = fftdim * fftdim; nucell = dm->nucell; nx = ny = nz = 5; write(1); char str[MAXLINE]; - if (count_words(fgets(str,MAXLINE,stdin)) >= 3){ + if (input == NULL) input = new UserInput(0); + input->read_stdin(str); + if (count_words(str) >= 3){ nx = atoi(strtok(str," \t\n\r\f")); ny = atoi(strtok(NULL," \t\n\r\f")); nz = atoi(strtok(NULL," \t\n\r\f")); @@ -36,7 +50,7 @@ Phonopy::Phonopy(DynMat *dynmat) memory->create(mass, nucell, "Phonopy:mass"); for (int i = 0; i < nucell; ++i){ - double m = 1./dm->M_inv_sqrt[i]; + double m = 1.0/dm->M_inv_sqrt[i]; mass[i] = m * m; } @@ -68,7 +82,7 @@ return; void Phonopy::write(int flag) { if (flag == 1){ // basic information - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + puts("================================================================================"); printf("Now to prepare the input files for phonopy.\n"); printf("The dimension of your present supercell is : %d x %d x %d.\n", dm->nx, dm->ny, dm->nz); printf("The size of the force constant matrix will be: %d x %d.\n", dm->npt*3, dm->npt*3); @@ -84,19 +98,18 @@ void Phonopy::write(int flag) } else if (flag == 4){ printf("Done!\nThe force constants information is extracted and written to FORCE_CONSTANTS,\n"); printf("the primitive cell is written to POSCAR.primitive, and the input file for\n"); - printf("phonopy band evaluation is written to band.conf.\n"); - printf("One should be able to obtain the phonon band structure after correcting\n"); - printf("the element names in POSCAR.primitive and band.conf by running\n"); - printf("`phonopy --readfc -c POSCAR.primitive -p band.conf`.\n"); - for (int ii = 0; ii < 80; ++ii) printf("-"); - printf("\n*** Remember to change the element names. ***\n"); -#ifdef UseSPG - for (int ii = 0; ii < 80; ++ii) printf("-"); -#endif + printf("phonopy band evaluation is written to band.conf.\n\n"); + printf("One should be able to obtain the phonon band structure after\n"); + printf(" 1) Correcting the `element names` in POSCAR.primitive and band.conf;\n"); + printf(" 2) Running `phonopy --readfc -c POSCAR.primitive -p band.conf`.\n\n"); + printf("Or the phonon density of states after\n"); + printf(" 1) Correcting the `element names` in POSCAR.primitive and mesh.conf;\n"); + printf(" 2) Running `phonopy --readfc -c POSCAR.primitive -p mesh.conf`.\n"); + puts("--------------------------------------------------------------------------------"); + printf("*** Remember to modify the `element names`. ***\n"); } else if (flag == 5){ - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); - + puts("================================================================================"); } return; } @@ -162,7 +175,9 @@ void Phonopy::phonopy() memory->destroy(out); // in POSCAR, atoms are sorted/aggregated by type, while for LAMMPS there is no such requirment - int type_id[nucell], num_type[nucell], ntype = 0; + int *type_id = new int[nucell]; + int *num_type = new int[nucell]; + int ntype = 0; for (int i = 0; i < nucell; ++i) num_type[i] = 0; for (int i = 0; i < nucell; ++i){ int ip = ntype; @@ -221,7 +236,14 @@ void Phonopy::phonopy() // write the primitive cell in POSCAR format fp = fopen("POSCAR.primitive", "w"); fprintf(fp, "Fix-phonon unit cell"); - for (int ip = 0; ip < ntype; ++ip) fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[ip]); + for (int ip = 0; ip < ntype; ++ip){ + for (int i = 0; i < nucell; ++i){ + if (dm->attyp[i] == type_id[ip]){ + fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]); + break; + } + } + } fprintf(fp, "\n1.\n"); int ndim = 0; for (int idim = 0; idim < 3; ++idim){ @@ -240,57 +262,48 @@ void Phonopy::phonopy() } fclose(fp); -#ifdef UseSPG - // Get high symmetry k-path - QNodes *q = new QNodes(); - kPath *kp = new kPath(dm, q); - kp->kpath(); -#endif - + // mesh.conf + fp = fopen("mesh.conf", "w"); + fprintf(fp, "# From Fix-phonon"); + for (int ip = 0; ip < ntype; ++ip){ + for (int i = 0; i < nucell; ++i){ + if (dm->attyp[i] == type_id[ip]){ + fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]); + break; + } + } + } + fprintf(fp, "\n\nATOM_NAME = "); + for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]); + fprintf(fp, "\nDIM = %d %d %d\n", nx, ny, nz); + fprintf(fp, "MP = 31 31 31\nFORCE_CONSTANTS = READ\n"); + fprintf(fp, "#FC_SYMMETRY = .TRUE.\n#SYMMETRY_TOLERANCE = 0.01\n"); + fclose(fp); + + // band.conf fp = fopen("band.conf", "w"); fprintf(fp, "# From Fix-phonon"); - for (int ip = 0; ip < ntype; ++ip) fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[ip]); + for (int ip = 0; ip < ntype; ++ip){ + for (int i = 0; i < nucell; ++i){ + if (dm->attyp[i] == type_id[ip]){ + fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]); + break; + } + } + } fprintf(fp, "\n\nATOM_NAME = "); for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]); - fprintf(fp, "\nDIM = %d %d %d\nBAND = ", nx, ny, nz); -#ifdef UseSPG - int nsect = q->qs.size(); - for (int i = 0; i < nsect; ++i){ - fprintf(fp, " %lg %lg %lg", q->qs[i][0], q->qs[i][1], q->qs[i][2]); - if (i+1 < nsect){ - double dq = 0.; - for (int j = 0; j < 3; ++j) dq += (q->qe[i][j] - q->qs[i+1][j]) * (q->qe[i][j] - q->qs[i+1][j]); - if (dq > ZERO) { - fprintf(fp, " %lg %lg %lg,", q->qe[i][0], q->qe[i][1], q->qe[i][2]); - } - } else if (i+1 == nsect){ - fprintf(fp, " %lg %lg %lg\n", q->qe[i][0], q->qe[i][1], q->qe[i][2]); - } - } -#endif - fprintf(fp, "\nBAND_POINTS = 21\nBAND_LABELS ="); -#ifdef UseSPG - for (int i = 0; i < q->ndstr.size(); ++i){ - std::size_t found = q->ndstr[i].find("{/Symbol G}"); - if (found != std::string::npos) q->ndstr[i].replace(found, found+11, "$\\Gamma$"); - found = q->ndstr[i].find("/"); - if (found != std::string::npos) q->ndstr[i].replace(found, found, " "); - fprintf(fp, " %s", q->ndstr[i].c_str()); - } -#endif - fprintf(fp, "\nFORCE_CONSTANTS = READ\nBAND_CONNECTION = .TRUE.\n"); + fprintf(fp, "\nDIM = %d %d %d\nBAND = AUTO\n", nx, ny, nz); + fprintf(fp, "BAND_POINTS = 21\nFORCE_CONSTANTS = READ\nBAND_CONNECTION = .TRUE.\n"); + fprintf(fp, "#FC_SYMMETRY = .TRUE.\n#SYMMETRY_TOLERANCE = 0.01\n"); // output info write(4); -#ifdef UseSPG - kp->show_path(); - delete kp; - delete q; -#endif write(5); - -return; + delete[] type_id; + delete[] num_type; + return; } /*------------------------------------------------------------------------------ @@ -304,7 +317,7 @@ int Phonopy::count_words(const char *line) strcpy(copy,line); char *ptr; - if (ptr = strchr(copy,'#')) *ptr = '\0'; + if ((ptr = strchr(copy,'#'))) *ptr = '\0'; if (strtok(copy," \t\n\r\f") == NULL) { memory->destroy(copy); @@ -314,7 +327,7 @@ int Phonopy::count_words(const char *line) while (strtok(NULL," \t\n\r\f")) n++; memory->destroy(copy); -return n; + return n; } /*----------------------------------------------------------------------------*/ #endif diff --git a/tools/phonon/phonopy.h b/tools/phonon/phonopy.h index 2f3006e791..e1200b8242 100644 --- a/tools/phonon/phonopy.h +++ b/tools/phonon/phonopy.h @@ -4,22 +4,16 @@ #ifndef PHONOPY_H #define PHONOPY_H -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "memory.h" -#include "qnodes.h" -#include "dynmat.h" -#include "global.h" +#include "zheevd.h" class Phonopy { public: - Phonopy(DynMat *); + Phonopy(class DynMat *); ~Phonopy(); private: - Memory *memory; - char str[MAXLINE]; + class UserInput *input; + class Memory *memory; int npt, fftdim2; // local variables int nx, ny, nz, nucell; // local variables int sysdim, fftdim; // local variables @@ -27,7 +21,7 @@ private: doublecomplex **FC_all; - DynMat *dm; + class DynMat *dm; void write(int); void get_my_FC(); void phonopy(); diff --git a/tools/phonon/timer.cpp b/tools/phonon/timer.cpp index 953257340d..b94bed40b1 100644 --- a/tools/phonon/timer.cpp +++ b/tools/phonon/timer.cpp @@ -1,4 +1,7 @@ #include "timer.h" + +#include + /* ----------------------------------------------------------------------------- * Initialization of time * -------------------------------------------------------------------------- */ diff --git a/tools/phonon/timer.h b/tools/phonon/timer.h index cd04dea56d..993a33f9f3 100644 --- a/tools/phonon/timer.h +++ b/tools/phonon/timer.h @@ -1,9 +1,7 @@ #ifndef TIMER_H #define TIMER_H -#include "stdio.h" -#include "stdlib.h" -#include "time.h" +#include class Timer { public: diff --git a/tools/phonon/tricubic/CMakeLists.txt b/tools/phonon/tricubic/CMakeLists.txt new file mode 100644 index 0000000000..d70b44c2b2 --- /dev/null +++ b/tools/phonon/tricubic/CMakeLists.txt @@ -0,0 +1,18 @@ + +# Support Linux from Ubuntu 20.04LTS onward, CentOS 7.x (with EPEL), +# macOS, MSVC 2019 (=Version 16) +cmake_minimum_required(VERSION 3.16) + +# set up project +project(tricubic VERSION 1.1 DESCRIPTION "Tricubic library" LANGUAGES CXX) +set(CMAKE_POSITION_INDEPENDENT_CODE ON) + +# hacks for MSVC to prevent lots of pointless warnings about "unsafe" functions +if(MSVC) + add_compile_options(/wd4244) + add_compile_options(/wd4267) + add_compile_definitions(_CRT_SECURE_NO_WARNINGS) +endif() + +add_library(tricubic STATIC tricubic.cpp) +target_include_directories(tricubic PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/tools/phonon/tricubic/LICENSE b/tools/phonon/tricubic/LICENSE new file mode 100644 index 0000000000..d60c31a97a --- /dev/null +++ b/tools/phonon/tricubic/LICENSE @@ -0,0 +1,340 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License. diff --git a/tools/phonon/tricubic/README.md b/tools/phonon/tricubic/README.md new file mode 100644 index 0000000000..0e8058339e --- /dev/null +++ b/tools/phonon/tricubic/README.md @@ -0,0 +1,28 @@ +# libtricubic + +This folder contains a slightly refactored version of version 1.0 of +the _libtricubic_ library, developed by François Lekien in 2004. + +The following paper explaining the method was published in 2005: + +> Lekien, F., & Marsden, J. (2005). _Tricubic interpolation in three +> dimensions._ International Journal for Numerical Methods in Engineering, +> 63(3), 455–471. doi:10.1002/nme.1296 + +Some additional notes and the full matrix can be found in the technical notes: + +> Lekien, F., Coulliette, C., & Marsden, J. (2004). _Tricubic Engine - Technical +> Notes and Full Matrix. + +The main article refers to the author's website +(http://gyre.cds.caltech.edu/pub/software/tricubic/) to download the code and +documentation. Unfortunately, this website no longer exists. Even the +archive.org snapshot is useless; [A single snapshot from December 3rd 2009 is +available](https://web.archive.org/web/20091203115835/http://gyre.cds.caltech.edu/pub/software/tricubic) +which seems like an FTP listing. No files are accessible. + +The source code was obtained from https://github.com/nbigaouette/libtricubic/ +Only the sources for the library were retained. No functional changes were +made, but some common programming conventions were applied, source files +merged, and build support for CMake added. + diff --git a/tools/phonon/tricubic/tricubic.cpp b/tools/phonon/tricubic/tricubic.cpp new file mode 100644 index 0000000000..974ed6aa52 --- /dev/null +++ b/tools/phonon/tricubic/tricubic.cpp @@ -0,0 +1,185 @@ + +#include "tricubic.h" + +#include + +static const int A[64][64] = { + { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 9,-9,-9, 9, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {-6, 6, 6,-6, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {-6, 6, 6,-6, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 4,-4,-4, 4, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}, + {-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 9,-9, 0, 0,-9, 9, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {-6, 6, 0, 0, 6,-6, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9, 0, 0,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0}, + { 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0}, + {-27,27,27,-27,27,-27,-27,27,-18,-9,18, 9,18, 9,-18,-9,-18,18,-9, 9,18,-18, 9,-9,-18,18,18,-18,-9, 9, 9,-9,-12,-6,-6,-3,12, 6, 6, 3,-12,-6,12, 6,-6,-3, 6, 3,-12,12,-6, 6,-6, 6,-3, 3,-8,-4,-4,-2,-4,-2,-2,-1}, + {18,-18,-18,18,-18,18,18,-18, 9, 9,-9,-9,-9,-9, 9, 9,12,-12, 6,-6,-12,12,-6, 6,12,-12,-12,12, 6,-6,-6, 6, 6, 6, 3, 3,-6,-6,-3,-3, 6, 6,-6,-6, 3, 3,-3,-3, 8,-8, 4,-4, 4,-4, 2,-2, 4, 4, 2, 2, 2, 2, 1, 1}, + {-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0}, + {18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6, 9,-9, 9,-9,-9, 9,-9, 9,12,-12,-12,12, 6,-6,-6, 6, 6, 3, 6, 3,-6,-3,-6,-3, 8, 4,-8,-4, 4, 2,-4,-2, 6,-6, 6,-6, 3,-3, 3,-3, 4, 2, 4, 2, 2, 1, 2, 1}, + {-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-6, 6,-6, 6, 6,-6, 6,-6,-8, 8, 8,-8,-4, 4, 4,-4,-3,-3,-3,-3, 3, 3, 3, 3,-4,-4, 4, 4,-2,-2, 2, 2,-4, 4,-4, 4,-2, 2,-2, 2,-2,-2,-2,-2,-1,-1,-1,-1}, + { 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {-6, 6, 0, 0, 6,-6, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 4,-4, 0, 0,-4, 4, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0}, + {-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0}, + {18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6,12,-12, 6,-6,-12,12,-6, 6, 9,-9,-9, 9, 9,-9,-9, 9, 8, 4, 4, 2,-8,-4,-4,-2, 6, 3,-6,-3, 6, 3,-6,-3, 6,-6, 3,-3, 6,-6, 3,-3, 4, 2, 2, 1, 4, 2, 2, 1}, + {-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-8, 8,-4, 4, 8,-8, 4,-4,-6, 6, 6,-6,-6, 6, 6,-6,-4,-4,-2,-2, 4, 4, 2, 2,-3,-3, 3, 3,-3,-3, 3, 3,-4, 4,-2, 2,-4, 4,-2, 2,-2,-2,-1,-1,-2,-2,-1,-1}, + { 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0}, + {-12,12,12,-12,12,-12,-12,12,-8,-4, 8, 4, 8, 4,-8,-4,-6, 6,-6, 6, 6,-6, 6,-6,-6, 6, 6,-6,-6, 6, 6,-6,-4,-2,-4,-2, 4, 2, 4, 2,-4,-2, 4, 2,-4,-2, 4, 2,-3, 3,-3, 3,-3, 3,-3, 3,-2,-1,-2,-1,-2,-1,-2,-1}, + { 8,-8,-8, 8,-8, 8, 8,-8, 4, 4,-4,-4,-4,-4, 4, 4, 4,-4, 4,-4,-4, 4,-4, 4, 4,-4,-4, 4, 4,-4,-4, 4, 2, 2, 2, 2,-2,-2,-2,-2, 2, 2,-2,-2, 2, 2,-2,-2, 2,-2, 2,-2, 2,-2, 2,-2, 1, 1, 1, 1, 1, 1, 1, 1} +}; + +static const char tricubic_version_stored[] = "1.1"; + +static int ijk2n(int i, int j, int k) { + return(i+4*j+16*k); +} + +static void point2xyz(int p, int *x, int *y, int *z) { + switch (p) { + case 0: *x=0; *y=0; *z=0; break; + case 1: *x=1; *y=0; *z=0; break; + case 2: *x=0; *y=1; *z=0; break; + case 3: *x=1; *y=1; *z=0; break; + case 4: *x=0; *y=0; *z=1; break; + case 5: *x=1; *y=0; *z=1; break; + case 6: *x=0; *y=1; *z=1; break; + case 7: *x=1; *y=1; *z=1; break; + default:*x=0; *y=0; *z=0; + } +} + +const char *tricubic_version(void) { + return(tricubic_version_stored); +} + +void tricubic_pointID2xyz(int id, int *x, int *y, int *z) { + point2xyz(id,x,y,z); +} + +void tricubic_pointID2xyz(int id, double *x, double *y, double *z) { + int x2,y2,z2; + point2xyz(id,&x2,&y2,&z2); + *x=(double)(x2); + *y=(double)(y2); + *z=(double)(z2); +} + +void tricubic_get_coeff_stacked(double a[64], double x[64]) { + int i,j; + for (i=0;i<64;i++) { + a[i]=(double)(0.0); + for (j=0;j<64;j++) { + a[i]+=A[i][j]*x[j]; + } + } +} + +void tricubic_get_coeff(double a[64], double f[8], double dfdx[8], double dfdy[8], double dfdz[8], double d2fdxdy[8], double d2fdxdz[8], double d2fdydz[8], double d3fdxdydz[8]) { + int i; + double x[64]; + for (i=0;i<8;i++) { + x[0+i]=f[i]; + x[8+i]=dfdx[i]; + x[16+i]=dfdy[i]; + x[24+i]=dfdz[i]; + x[32+i]=d2fdxdy[i]; + x[40+i]=d2fdxdz[i]; + x[48+i]=d2fdydz[i]; + x[56+i]=d3fdxdydz[i]; + } + tricubic_get_coeff_stacked(a,x); +} + +double tricubic_eval(double a[64], double x, double y, double z) { + int i,j,k; + double ret=(double)(0.0); + /* TRICUBIC EVAL + This is the short version of tricubic_eval. It is used to compute + the value of the function at a given point (x,y,z). To compute + partial derivatives of f, use the full version with the extra args. + */ + for (i=0;i<4;i++) { + for (j=0;j<4;j++) { + for (k=0;k<4;k++) { + ret+=a[ijk2n(i,j,k)]*pow(x,i)*pow(y,j)*pow(z,k); + } + } + } + return(ret); +} + +double tricubic_eval(double a[64], double x, double y, double z, int derx, int dery, int derz) { + int i,j,k; + double ret=(double)(0.0); + double cont; + int w; + /* TRICUBIC_EVAL + The full version takes 3 extra integers args that allows to evaluate + any partial derivative of f at the point + derx=dery=derz=0 => f + derx=2 dery=derz=0 => d2f/dx2 + derx=dery=derz=1 =? d3f/dxdydz + NOTICE that (derx>3)||(dery>3)||(derz>3) => returns 0.0 + this computes \frac{\partial ^{derx+dery+derz} d}{\partial x ^{derx} \partial y ^{dery} \partial z ^{derz}} + */ + for (i=derx;i<4;i++) { + for (j=dery;j<4;j++) { + for (k=derz;k<4;k++) { + cont=a[ijk2n(i,j,k)]*pow(x,i-derx)*pow(y,j-dery)*pow(z,k-derz); + for (w=0;w Date: Sat, 29 Apr 2023 03:55:54 -0400 Subject: [PATCH 21/41] update so it can be built either standalone or as subdir with LAMMPS --- tools/phonon/CMakeLists.txt | 21 +++++++++++++++++++++ tools/phonon/version.h.in | 1 + 2 files changed, 22 insertions(+) create mode 100644 tools/phonon/version.h.in diff --git a/tools/phonon/CMakeLists.txt b/tools/phonon/CMakeLists.txt index 283b995590..02ccee25a3 100644 --- a/tools/phonon/CMakeLists.txt +++ b/tools/phonon/CMakeLists.txt @@ -47,6 +47,10 @@ add_executable(phana ) target_include_directories(phana PUBLIC $) +if(NOT LAMMPS_DIR) + set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../../cmake/Modules) + set(LAMMPS_THIRDPARTY_URL "https://download.lammps.org/thirdparty") +endif() find_package(FFTW3) if(FFTW3_FOUND) target_compile_definitions(phana PRIVATE FFTW3) @@ -56,12 +60,29 @@ endif() # build bundeled libraries add_subdirectory(tricubic) +# standalone build must build our own version of linalg +if(NOT LAMMPS_DIR) + if(NOT USE_INTERNAL_LINALG) + find_package(LAPACK) + find_package(BLAS) + endif() + if(NOT LAPACK_FOUND OR NOT BLAS_FOUND OR USE_INTERNAL_LINALG) + file(GLOB LINALG_SOURCES CONFIGURE_DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/../../lib/linalg/[^.]*.cpp) + add_library(linalg STATIC ${LINALG_SOURCES}) + set(BLAS_LIBRARIES "$") + set(LAPACK_LIBRARIES "$") + else() + list(APPEND LAPACK_LIBRARIES ${BLAS_LIBRARIES}) + endif() +endif() + option(USE_SPGLIB "Download and use spglib for phonon DOS and other optional properties" ON) if(USE_SPGLIB) set(SPGLIB_URL "https://github.com/spglib/spglib/archive/refs/tags/v1.11.2.1.tar.gz" CACHE STRING "URL for spglib v1.x tarball") set(SPGLIB_MD5 "3089782bc85b5034dd4765a18ee70bc7" CACHE STRING "MD5 checksum for spglib tarball") mark_as_advanced(SPGLIB_URL) mark_as_advanced(SPGLIB_MD5) + include(LAMMPSUtils) GetFallbackURL(SPGLIB_URL SPGLIB_FALLBACK) include(ExternalProject) diff --git a/tools/phonon/version.h.in b/tools/phonon/version.h.in new file mode 100644 index 0000000000..1f6fb01373 --- /dev/null +++ b/tools/phonon/version.h.in @@ -0,0 +1 @@ +#define VERSION @CMAKE_PROJECT_VERSION@ From a2b74c846648a6b956aa66915a79340f1f6f2ee2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 20:23:13 -0400 Subject: [PATCH 22/41] small doc update --- doc/src/Tools.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst index 2eae1aabc7..bccb9abe73 100644 --- a/doc/src/Tools.rst +++ b/doc/src/Tools.rst @@ -883,9 +883,9 @@ dependencies and redirects the download to the local cache. phonon tool ------------------------ -The phonon subdirectory contains a post-processing tool useful for -analyzing the output of the :doc:`fix phonon ` command in -the PHONON package. +The phonon subdirectory contains a post-processing tool, *phana*, useful +for analyzing the output of the :doc:`fix phonon ` command +in the PHONON package. See the README file for instruction on building the tool and what library it needs. And see the examples/PACKAGES/phonon directory From 73b8bb8617cdddc2a2a812090d40bb26e2a9b539 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 20:41:44 -0400 Subject: [PATCH 23/41] fix bug that breaks compilation with -DLAMMPS_BIGBIG --- src/compute_count_type.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index feadcb04ba..4fd814cdff 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -144,7 +144,7 @@ void ComputeCountType::compute_vector() // bond types can be negative for SHAKE else if (mode == BOND) { - int **bond_atom = atom->bond_atom; + tagint **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; int *num_bond = atom->num_bond; int *mask = atom->mask; From b38544a9e86612e130690dcd56460c3d50802a0a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 20:43:05 -0400 Subject: [PATCH 24/41] apply clang-format --- src/compute_count_type.cpp | 50 ++++++++++++++++++++------------------ src/compute_count_type.h | 2 +- 2 files changed, 27 insertions(+), 25 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 4fd814cdff..786506280e 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -23,7 +23,7 @@ using namespace LAMMPS_NS; -enum{ATOM,BOND}; +enum { ATOM, BOND }; /* ---------------------------------------------------------------------- */ @@ -33,9 +33,12 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( // process args - if (strcmp(arg[3],"atom") == 0) mode = ATOM; - else if (strcmp(arg[3],"bond") == 0) mode = BOND; - else error->all(FLERR, "Invalid compute count/type keyword {}",arg[3]); + if (strcmp(arg[3], "atom") == 0) + mode = ATOM; + else if (strcmp(arg[3], "bond") == 0) + mode = BOND; + else + error->all(FLERR, "Invalid compute count/type keyword {}", arg[3]); if (mode == ATOM) { vector_flag = 1; @@ -49,10 +52,10 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( } if (mode == BOND && !atom->avec->bonds_allow) - error->all(FLERR,"Cannot use compute count/type bond command with no bonds allowed"); + error->all(FLERR, "Cannot use compute count/type bond command with no bonds allowed"); // output vector - + vector = new double[size_vector]; // work vectors @@ -86,13 +89,13 @@ double ComputeCountType::compute_scalar() int **bond_type = atom->bond_type; int nlocal = atom->nlocal; - int m,nbond; + int m, nbond; int count = 0; // count broken bonds with bond_type = 0 // ignore group setting since 2 atoms in a broken bond // can be arbitrarily far apart - + for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; for (m = 0; m < nbond; m++) @@ -101,14 +104,13 @@ double ComputeCountType::compute_scalar() // sum across procs as bigint, then convert to double // correct for double counting if newton_bond off - + bigint bcount_me = count; bigint bcount; MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); if (force->newton_bond == 0) bcount /= 2; - if (bcount > MAXDOUBLEINT) - error->all(FLERR,"Compute count/type overflow"); + if (bcount > MAXDOUBLEINT) error->all(FLERR, "Compute count/type overflow"); scalar = bcount; return scalar; } @@ -132,12 +134,11 @@ void ComputeCountType::compute_vector() for (int m = 0; m < ntypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - count[type[i]-1]++; + if (mask[i] & groupbit) count[type[i] - 1]++; nvec = ntypes; } - + // count bonds by type // both atoms in bond must be in group to be counted // skip type = 0 bonds, they are counted by compute_scalar @@ -150,11 +151,11 @@ void ComputeCountType::compute_vector() int *mask = atom->mask; int nlocal = atom->nlocal; int nbondtypes = atom->nbondtypes; - - int j,m,nbond,itype; + + int j, m, nbond, itype; int flag = 0; for (int m = 0; m < nbondtypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; for (m = 0; m < nbond; m++) { @@ -166,17 +167,19 @@ void ComputeCountType::compute_vector() flag = 1; continue; } - + if ((mask[i] & groupbit) && (mask[j] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); + if (flagany) error->all(FLERR, "Missing bond atom in compute count/type"); nvec = nbondtypes; } @@ -186,11 +189,10 @@ void ComputeCountType::compute_vector() for (int m = 0; m < nvec; m++) bcount_me[m] = count[m]; MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world); - if (force->newton_bond == 0) + if (force->newton_bond == 0) for (int m = 0; m < nvec; m++) bcount[m] /= 2; for (int m = 0; m < nvec; m++) - if (bcount[m] > MAXDOUBLEINT) - error->all(FLERR,"Compute count/type overflow"); + if (bcount[m] > MAXDOUBLEINT) error->all(FLERR, "Compute count/type overflow"); for (int m = 0; m < nvec; m++) vector[m] = bcount[m]; } diff --git a/src/compute_count_type.h b/src/compute_count_type.h index 0ed7a22467..27b1669b1d 100644 --- a/src/compute_count_type.h +++ b/src/compute_count_type.h @@ -34,7 +34,7 @@ class ComputeCountType : public Compute { protected: int mode; - + int *count; bigint *bcount_me; bigint *bcount; From ff29ef7d3160c9597b9fb637619982afcf8e0adc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 20:57:32 -0400 Subject: [PATCH 25/41] use more obvious representation of 2^53 --- src/lmptype.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lmptype.h b/src/lmptype.h index 891d5bdf89..769f730d6a 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -95,8 +95,8 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT_MAX #define MAXBIGINT INT64_MAX -#define MAXDOUBLEINT 9007199254740992 // 2^53 - +#define MAXDOUBLEINT (1L<<53) + #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT #define MPI_LMP_BIGINT MPI_LL @@ -133,7 +133,7 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT64_MAX #define MAXBIGINT INT64_MAX -#define MAXDOUBLEINT 9007199254740992 // 2^53 +#define MAXDOUBLEINT (1L<<53) #define MPI_LMP_TAGINT MPI_LL #define MPI_LMP_IMAGEINT MPI_LL From c8d5b9e4d072c46c1795e696e26fc09b48016510 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 23:52:04 -0400 Subject: [PATCH 26/41] avoid uninitialized pointers --- src/compute_count_type.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 786506280e..0cd109a1eb 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -27,9 +27,10 @@ enum { ATOM, BOND }; /* ---------------------------------------------------------------------- */ -ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) +ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), count(nullptr), bcount(nullptr), bcount_me(nullptr) { - if (narg != 4) error->all(FLERR, "Illegal compute count/type command"); + if (narg != 4) error->all(FLERR, "Incorrect number of args for compute count/type command"); // process args @@ -51,7 +52,7 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( extvector = 1; } - if (mode == BOND && !atom->avec->bonds_allow) + if ((mode == BOND) && !atom->avec->bonds_allow) error->all(FLERR, "Cannot use compute count/type bond command with no bonds allowed"); // output vector From 12135bac770ecc591b05a31ad437a72038553b14 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 10:22:37 -0400 Subject: [PATCH 27/41] add unit tests for compute count/type --- unittest/commands/test_compute_global.cpp | 54 +++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index d7573d8d16..4efc3f1b91 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -294,6 +294,60 @@ TEST_F(ComputeGlobalTest, Reduction) EXPECT_DOUBLE_EQ(rep[2], 26); EXPECT_DOUBLE_EQ(rep[3], max[0]); } + +TEST_F(ComputeGlobalTest, Counts) +{ + if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP(); + + BEGIN_HIDE_OUTPUT(); + command("pair_style lj/cut 10.0"); + command("pair_coeff * * 0.01 3.0"); + command("bond_style harmonic"); + command("bond_coeff * 100.0 1.5"); + + command("variable t1 atom type==1"); + command("variable t2 atom type==2"); + command("variable t3 atom type==3"); + command("variable t4 atom type==4"); + command("variable t5 atom type==5"); + command("compute asum all reduce sum v_t1 v_t2 v_t3 v_t4 v_t5"); + command("compute acnt all count/type atom"); + command("compute bcnt all count/type bond"); + command("thermo_style custom c_asum[*] c_acnt[*]"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + + auto asum = get_vector("asum"); + auto acnt = get_vector("acnt"); + auto bcnt = get_vector("bcnt"); + auto bbrk = get_scalar("bcnt"); + + EXPECT_DOUBLE_EQ(bbrk, 0.0); + + EXPECT_DOUBLE_EQ(asum[0], acnt[0]); + EXPECT_DOUBLE_EQ(asum[1], acnt[1]); + EXPECT_DOUBLE_EQ(asum[2], acnt[2]); + EXPECT_DOUBLE_EQ(asum[3], acnt[3]); + EXPECT_DOUBLE_EQ(asum[4], acnt[4]); + + EXPECT_DOUBLE_EQ(bcnt[0], 3.0); + EXPECT_DOUBLE_EQ(bcnt[1], 6.0); + EXPECT_DOUBLE_EQ(bcnt[2], 3.0); + EXPECT_DOUBLE_EQ(bcnt[3], 2.0); + EXPECT_DOUBLE_EQ(bcnt[4], 10.0); + + BEGIN_HIDE_OUTPUT(); + command("delete_bonds all bond 3"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + + EXPECT_DOUBLE_EQ(bbrk, 0.0); // should be 3 + EXPECT_DOUBLE_EQ(bcnt[0], 3.0); + EXPECT_DOUBLE_EQ(bcnt[1], 6.0); + EXPECT_DOUBLE_EQ(bcnt[2], 3.0); // should be 0 + EXPECT_DOUBLE_EQ(bcnt[3], 2.0); + EXPECT_DOUBLE_EQ(bcnt[4], 10.0); +} } // namespace LAMMPS_NS int main(int argc, char **argv) From eacb420e211452392442b34b26add5693dd3c092 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 10:34:36 -0400 Subject: [PATCH 28/41] whitespace --- doc/src/compute_count_type.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index 9ab0b1009f..8cbd6aec9b 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -13,7 +13,7 @@ Syntax * ID, group-ID are documented in :doc:`compute ` command * count/type = style name of this compute command * mode = {atom} or {bond} - + Examples """""""" From d2bb1b420da791f73b141fc1ad52027b17846d16 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 11:49:13 -0400 Subject: [PATCH 29/41] simplify and cleanup --- unittest/commands/test_compute_global.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index 4efc3f1b91..7e5b3d4483 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -300,10 +300,8 @@ TEST_F(ComputeGlobalTest, Counts) if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP(); BEGIN_HIDE_OUTPUT(); - command("pair_style lj/cut 10.0"); - command("pair_coeff * * 0.01 3.0"); - command("bond_style harmonic"); - command("bond_coeff * 100.0 1.5"); + command("pair_style zero 10.0"); + command("pair_coeff * *"); command("variable t1 atom type==1"); command("variable t2 atom type==2"); @@ -341,6 +339,8 @@ TEST_F(ComputeGlobalTest, Counts) command("run 0 post no"); END_HIDE_OUTPUT(); + bcnt = get_vector("bcnt"); + bbrk = get_scalar("bcnt"); EXPECT_DOUBLE_EQ(bbrk, 0.0); // should be 3 EXPECT_DOUBLE_EQ(bcnt[0], 3.0); EXPECT_DOUBLE_EQ(bcnt[1], 6.0); From aa4447413ae32c3c51db83cb75facb8eafcb0d12 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 2 May 2023 13:39:37 -0600 Subject: [PATCH 30/41] expand to angles, dihedrals, impropers --- doc/src/Howto_broken_bonds.rst | 35 ++-- doc/src/compute_count_type.rst | 68 +++++-- src/compute_count_type.cpp | 330 ++++++++++++++++++++++++++------- src/compute_count_type.h | 6 + 4 files changed, 335 insertions(+), 104 deletions(-) diff --git a/doc/src/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index 4fd0bcc30f..282d2b32a4 100644 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -1,9 +1,9 @@ Broken Bonds ============ -Typically, bond interactions persist for the duration of a simulation -in LAMMPS. However, some commands allow bonds to break, including the -following: +Typically, molecular bond interactions persist for the duration of a +simulation in LAMMPS. However, some commands break bonds dynamically, +including the following: * :doc:`bond_style quartic ` * :doc:`fix bond/break ` @@ -14,39 +14,40 @@ A bond can break if it is stretched beyond a user-defined threshold or more generally if other criteria are met. For the quartic bond style, when a bond is broken its bond type is set -to 0 and pairwise forces between the two atoms in the broken bond are -"turned on". Angles, dihedrals, etc cannot be defined for the system -when :doc:`bond_style quartic ` is used. +to 0 to effectively break it and pairwise forces between the two atoms +in the broken bond are "turned on". Angles, dihedrals, etc cannot be +defined for a system when :doc:`bond_style quartic ` is +used. The :doc:`fix bond/break ` and :doc:`fix bond/react ` commands allow breaking of bonds within a molecular -topology with also defines angles, dihedrals, etc. These fixes will -update internal topology data structures when bonds are broken, so -that the appropriate angle, dihederal, etc interactions are also -turned off. They will also trigger a rebuild of the neighbor list +topology with may also define angles, dihedrals, etc. These commands +update internal topology data structures to remove broken bonds, as +well as the appropriate angle, dihederal, etc interactions which +include the bond. They also trigger a rebuild of the neighbor list when this occurs, to turn on the appropriate pairwise forces. In the BPM package, one can either turn off all pair interactions between bonded particles or leave them on, overlaying pair forces on top of bond forces. To remove pair forces, the special bond list is -dynamically updated. More details can be found on the :doc:`Howto BPM +dynamically updated. More details can be found on the :doc:`Howto BPM ` page. Note that when bonds are dumped to a file via the :doc:`dump local ` command, bonds with type 0 are not included. -The :doc:`delete_bonds ` command can also be used to -query the status of broken bonds (type = 0) or permanently delete -them, e.g.: +The :doc:`delete_bonds ` command can be used to query +the status of broken bonds with type = 0 or permanently delete them, +e.g.: .. code-block:: LAMMPS delete_bonds all stats delete_bonds all bond 0 remove -The compute :doc:`count/type bond ` command -tallies the current number of bonds for each bond type. It also -tallies broken bonds with type = 0. +The compute :doc:`count/type ` command tallies the +current number of bonds (or angles, etc) for each bond (angle, etc) +type. It also tallies broken bonds with type = 0. The compute :doc:`nbond/atom ` command tallies the current number of bonds each atom is part of, excluding broken bonds diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index 9ab0b1009f..e9f531ef64 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -12,7 +12,7 @@ Syntax * ID, group-ID are documented in :doc:`compute ` command * count/type = style name of this compute command -* mode = {atom} or {bond} +* mode = {atom} or {bond} or {angle} or {dihedral} or {improper} Examples """""""" @@ -26,28 +26,46 @@ Description """"""""""" Define a computation that counts the current number of atoms for each -atom type or the number of bonds for each bond type. The former can -be useful if atoms are added to or deleted from the system in random -ways, e.g. via the :doc:`fix deposit `, :doc:`fix pour -`, or :doc:`fix evaporate ` commands. The -latter can be useful in reactive simulations where molecular bonds are -broken or created. +atom type. Or the number of bonds (angles, dihedrals, impropers) for +each bond (angle, dihedral, improper) type. -Note that for this command, bonds are defined as the topological kind -enumerated in a data file, initially read by the :doc:`read_data +The former can be useful if atoms are added to or deleted from the +system in random ways, e.g. via the :doc:`fix deposit `, +:doc:`fix pour `, or :doc:`fix evaporate ` +commands. The latter can be useful in reactive simulations where +molecular bonds are broken or created, as well as angles, dihedrals, +impropers. + +Note that for this command, bonds (angles, etc) are the topological +kind enumerated in a data file, initially read by the :doc:`read_data ` command or defined by the :doc:`molecule ` command. They do not refer to implicit bonds defined on-the-fly by bond-order or reactive pair styles based on the current conformation of small clusters of atoms. -These commands can create and/or break topological bonds: +These commands can turn off topological bonds (angles, etc) by setting +their bond (angle, etc) types to negative values. This command +includes the turned-off bonds (angles, etc) in the count for each +type: +* :doc:`fix shake ` +* :doc:`delete_bonds ` + +These commands can create and/or break topological bonds (angles, +etc). In the case of breaking, they remove the bond (angle, etc) from +the system, so that they no longer exist (:doc:`bond_style quartic +` is the exception, see the discussion below). Thus +they are not included in the counts for each type: + +* :doc:`delete_bonds remove ` +* :doc:`bond_style quartic ` * :doc:`fix bond/react ` * :doc:`fix bond/create ` * :doc:`fix bond/break ` -* :doc:`bond_style quartic ` * :doc:`BPM package ` bond styles +---------- + If the {mode} setting is {atom} then the count of atoms for each atom type is tallied. Only atoms in the specified group are counted. @@ -56,11 +74,22 @@ type is tallied. Only bonds with both atoms in the specified group are counted. For {mode} = {bond}, broken bonds with a bond type of zero are also -counted. Some of the commands listed above which break bonds, do this -by setting their types to zero. See the :doc:`Howto broken bonds -` doc page for details. Note that the group -setting is ignored for broken bonds; all broken bonds in the system -are counted. +counted. The :doc:`bond_style quartic ` breaks a bond +by doing this. See the :doc:`Howto broken bonds ` +doc page for more details. Note that the group setting is ignored for +broken bonds; all broken bonds in the system are counted. + +If the {mode} setting is {angle} then the count of angles for each +angle type is tallied. Only angles with all 3 atoms in the specified +group are counted. + +If the {mode} setting is {dihedral} then the count of dihedrals for +each dihedral type is tallied. Only dihedrals with all 4 atoms in the +specified group are counted. + +If the {mode} setting is {improper} then the count of impropers for +each improper type is tallied. Only impropers with all 4 atoms in the +specified group are counted. ---------- @@ -68,11 +97,12 @@ Output info """"""""""" This compute calculates a global vector of counts. If the mode is -{atom}, the vector length is the number of atom types. If the mode is -{bond}, the vector length is the number of bond types. +{atom} or {bond} or {angle} or {dihedral} or {improper}, then the +vector length is the number of atom types or bond types or angle types +or dihedral types or improper types, respectively. If the mode is {bond} this compute also calculates a global scalar -which is the number of broken bonds. +which is the number of broken bonds with type = 0, as explained above. These values can be used by any command that uses global scalar or vector values from a compute as input. See the :doc:`Howto output diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index feadcb04ba..36aa21da70 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -14,7 +14,6 @@ #include "compute_count_type.h" #include "atom.h" -#include "atom_vec.h" #include "domain.h" #include "error.h" #include "force.h" @@ -23,7 +22,7 @@ using namespace LAMMPS_NS; -enum{ATOM,BOND}; +enum{ATOM,BOND,ANGLE,DIHEDRAL,IMPROPER}; /* ---------------------------------------------------------------------- */ @@ -35,8 +34,24 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( if (strcmp(arg[3],"atom") == 0) mode = ATOM; else if (strcmp(arg[3],"bond") == 0) mode = BOND; + else if (strcmp(arg[3],"angle") == 0) mode = ANGLE; + else if (strcmp(arg[3],"dihedral") == 0) mode = DIHEDRAL; + else if (strcmp(arg[3],"improper") == 0) mode = IMPROPER; else error->all(FLERR, "Invalid compute count/type keyword {}",arg[3]); + // error check + + if (mode == BOND && !atom->nbondtypes) + error->all(FLERR,"Compute count/type bond command with no bonds defined"); + if (mode == ANGLE && !atom->nangletypes) + error->all(FLERR,"Compute count/type bond command with no angles defined"); + if (mode == DIHEDRAL && !atom->ndihedraltypes) + error->all(FLERR,"Compute count/type dihedral command with no dihedrals defined"); + if (mode == IMPROPER && !atom->nimpropertypes) + error->all(FLERR,"Compute count/type improper command with no impropers defined"); + + // set vector lengths + if (mode == ATOM) { vector_flag = 1; size_vector = atom->ntypes; @@ -46,11 +61,20 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( size_vector = atom->nbondtypes; extscalar = 1; extvector = 1; + } else if (mode == ANGLE) { + vector_flag = 1; + size_vector = atom->nangletypes; + extvector = 1; + } else if (mode == DIHEDRAL) { + vector_flag = 1; + size_vector = atom->ndihedraltypes; + extvector = 1; + } else if (mode == IMPROPER) { + vector_flag = 1; + size_vector = atom->nimpropertypes; + extvector = 1; } - if (mode == BOND && !atom->avec->bonds_allow) - error->all(FLERR,"Cannot use compute count/type bond command with no bonds allowed"); - // output vector vector = new double[size_vector]; @@ -121,76 +145,246 @@ void ComputeCountType::compute_vector() int nvec; - // count atoms by type - // atom must be in group to be counted - - if (mode == ATOM) { - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int ntypes = atom->ntypes; - - for (int m = 0; m < ntypes; m++) count[m] = 0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - count[type[i]-1]++; - - nvec = ntypes; - } + if (mode == ATOM) nvec = count_atoms(); + else if (mode == BOND) nvec = count_bonds(); + else if (mode == ANGLE) nvec = count_angles(); + else if (mode == DIHEDRAL) nvec = count_dihedrals(); + else if (mode == IMPROPER) nvec = count_impropers(); - // count bonds by type - // both atoms in bond must be in group to be counted - // skip type = 0 bonds, they are counted by compute_scalar - // bond types can be negative for SHAKE - - else if (mode == BOND) { - int **bond_atom = atom->bond_atom; - int **bond_type = atom->bond_type; - int *num_bond = atom->num_bond; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int nbondtypes = atom->nbondtypes; - - int j,m,nbond,itype; - int flag = 0; - for (int m = 0; m < nbondtypes; m++) count[m] = 0; - - for (int i = 0; i < nlocal; i++) { - nbond = num_bond[i]; - for (m = 0; m < nbond; m++) { - itype = bond_type[i][m]; - if (itype == 0) continue; - - j = atom->map(bond_atom[i][m]); - if (j < 0) { - flag = 1; - continue; - } - - if ((mask[i] & groupbit) && (mask[j] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; - } - } - } - - int flagany; - MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); - - nvec = nbondtypes; - } - // sum across procs as bigint, then convert to double - // correct for double counting if newton_bond off + // correct for multiple counting if newton_bond off for (int m = 0; m < nvec; m++) bcount_me[m] = count[m]; MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world); - if (force->newton_bond == 0) - for (int m = 0; m < nvec; m++) bcount[m] /= 2; - + + if (force->newton_bond == 0) { + if (mode == BOND) + for (int m = 0; m < nvec; m++) bcount[m] /= 2; + else if (mode == ANGLE) + for (int m = 0; m < nvec; m++) bcount[m] /= 3; + if (mode == DIHEDRAL || mode == IMPROPER) + for (int m = 0; m < nvec; m++) bcount[m] /= 4; + } + for (int m = 0; m < nvec; m++) if (bcount[m] > MAXDOUBLEINT) error->all(FLERR,"Compute count/type overflow"); for (int m = 0; m < nvec; m++) vector[m] = bcount[m]; } + +/* ---------------------------------------------------------------------- + count atoms by type + atom must be in group to be counted +---------------------------------------------------------------------- */ + +int ComputeCountType::count_atoms() +{ + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int ntypes = atom->ntypes; + + for (int m = 0; m < ntypes; m++) count[m] = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + count[type[i]-1]++; + + return ntypes; +} + +/* ---------------------------------------------------------------------- + count bonds by type + both atoms in bond must be in group to be counted + skip type = 0 bonds, they are counted by compute_scalar() + bond types can be negative, count them as if positive +---------------------------------------------------------------------- */ + +int ComputeCountType::count_bonds() +{ + int **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *num_bond = atom->num_bond; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nbondtypes = atom->nbondtypes; + + int j,m,nbond,itype; + int flag = 0; + for (int m = 0; m < nbondtypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + nbond = num_bond[i]; + for (m = 0; m < nbond; m++) { + itype = bond_type[i][m]; + if (itype == 0) continue; + + j = atom->map(bond_atom[i][m]); + if (j < 0) { + flag = 1; + continue; + } + + if ((mask[i] & groupbit) && (mask[j] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + } + + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); + + return nbondtypes; +} + +/* ---------------------------------------------------------------------- + count angles by type + all 3 atoms in angle must be in group to be counted + angle types can be negative, count them as if positive +---------------------------------------------------------------------- */ + +int ComputeCountType::count_angles() +{ + int **angle_atom1 = atom->angle_atom1; + int **angle_atom2 = atom->angle_atom2; + int **angle_atom3 = atom->angle_atom3; + int **angle_type = atom->angle_type; + int *num_angle = atom->num_angle; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nangletypes = atom->nangletypes; + + int j1,j2,j3,m,nangle,itype; + int flag = 0; + for (int m = 0; m < nangletypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + nangle = num_angle[i]; + for (m = 0; m < nangle; m++) { + itype = angle_type[i][m]; + + j1 = atom->map(angle_atom1[i][m]); + j2 = atom->map(angle_atom2[i][m]); + j3 = atom->map(angle_atom3[i][m]); + if (j1 < 0 || j2 < 0 || j3 < 0) { + flag = 1; + continue; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && + (mask[j3] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + } + + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing angle atom in compute count/type"); + + return nangletypes; +} + +/* ---------------------------------------------------------------------- + count dihedrals by type + all 4 atoms in dihedral must be in group to be counted + dihedral types can be negative, count them as if positive +---------------------------------------------------------------------- */ + +int ComputeCountType::count_dihedrals() +{ + int **dihedral_atom1 = atom->dihedral_atom1; + int **dihedral_atom2 = atom->dihedral_atom2; + int **dihedral_atom3 = atom->dihedral_atom3; + int **dihedral_atom4 = atom->dihedral_atom4; + int **dihedral_type = atom->dihedral_type; + int *num_dihedral = atom->num_dihedral; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int ndihedraltypes = atom->ndihedraltypes; + + int j1,j2,j3,j4,m,ndihedral,itype; + int flag = 0; + for (int m = 0; m < ndihedraltypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + ndihedral = num_dihedral[i]; + for (m = 0; m < ndihedral; m++) { + itype = dihedral_type[i][m]; + + j1 = atom->map(dihedral_atom1[i][m]); + j2 = atom->map(dihedral_atom2[i][m]); + j3 = atom->map(dihedral_atom3[i][m]); + j4 = atom->map(dihedral_atom4[i][m]); + if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) { + flag = 1; + continue; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && + (mask[j3] & groupbit) && (mask[j4] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + } + + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing dihedral atom in compute count/type"); + + return ndihedraltypes; +} + +/* ---------------------------------------------------------------------- + count impropers by type + all 4 atoms in improper must be in group to be counted + improper types can be negative, count them as if positive +---------------------------------------------------------------------- */ + +int ComputeCountType::count_impropers() +{ + int **improper_atom1 = atom->improper_atom1; + int **improper_atom2 = atom->improper_atom2; + int **improper_atom3 = atom->improper_atom3; + int **improper_atom4 = atom->improper_atom4; + int **improper_type = atom->improper_type; + int *num_improper = atom->num_improper; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nimpropertypes = atom->nimpropertypes; + + int j1,j2,j3,j4,m,nimproper,itype; + int flag = 0; + for (int m = 0; m < nimpropertypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + nimproper = num_improper[i]; + for (m = 0; m < nimproper; m++) { + itype = improper_type[i][m]; + + j1 = atom->map(improper_atom1[i][m]); + j2 = atom->map(improper_atom2[i][m]); + j3 = atom->map(improper_atom3[i][m]); + j4 = atom->map(improper_atom4[i][m]); + if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) { + flag = 1; + continue; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && + (mask[j3] & groupbit) && (mask[j4] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + } + + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing improper atom in compute count/type"); + + return nimpropertypes; +} diff --git a/src/compute_count_type.h b/src/compute_count_type.h index 0ed7a22467..68a51cb42a 100644 --- a/src/compute_count_type.h +++ b/src/compute_count_type.h @@ -38,6 +38,12 @@ class ComputeCountType : public Compute { int *count; bigint *bcount_me; bigint *bcount; + + int count_atoms(); + int count_bonds(); + int count_angles(); + int count_dihedrals(); + int count_impropers(); }; } // namespace LAMMPS_NS From 8a9091595df80f2834cb7452e8d1cc0a8cd1a7b0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 15:56:45 -0400 Subject: [PATCH 31/41] Revert "use more obvious representation of 2^53" This reverts commit ff29ef7d3160c9597b9fb637619982afcf8e0adc as it conflicts with how MSVC interprets 1L (it would require 1LL which creates issues with other compilers). --- src/lmptype.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lmptype.h b/src/lmptype.h index 769f730d6a..891d5bdf89 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -95,8 +95,8 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT_MAX #define MAXBIGINT INT64_MAX -#define MAXDOUBLEINT (1L<<53) - +#define MAXDOUBLEINT 9007199254740992 // 2^53 + #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT #define MPI_LMP_BIGINT MPI_LL @@ -133,7 +133,7 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT64_MAX #define MAXBIGINT INT64_MAX -#define MAXDOUBLEINT (1L<<53) +#define MAXDOUBLEINT 9007199254740992 // 2^53 #define MPI_LMP_TAGINT MPI_LL #define MPI_LMP_IMAGEINT MPI_LL From 3119434932f05b1276ece2c7d17bfbd1076f44fd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 15:59:34 -0400 Subject: [PATCH 32/41] apply clang-format --- src/compute_count_type.cpp | 147 +++++++++++++++++++++---------------- 1 file changed, 82 insertions(+), 65 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index ed0467c435..5c27dd39ca 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -22,7 +22,7 @@ using namespace LAMMPS_NS; -enum{ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER}; +enum { ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER }; /* ---------------------------------------------------------------------- */ @@ -33,23 +33,29 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : // process args - if (strcmp(arg[3],"atom") == 0) mode = ATOM; - else if (strcmp(arg[3],"bond") == 0) mode = BOND; - else if (strcmp(arg[3],"angle") == 0) mode = ANGLE; - else if (strcmp(arg[3],"dihedral") == 0) mode = DIHEDRAL; - else if (strcmp(arg[3],"improper") == 0) mode = IMPROPER; - else error->all(FLERR, "Invalid compute count/type keyword {}",arg[3]); + if (strcmp(arg[3], "atom") == 0) + mode = ATOM; + else if (strcmp(arg[3], "bond") == 0) + mode = BOND; + else if (strcmp(arg[3], "angle") == 0) + mode = ANGLE; + else if (strcmp(arg[3], "dihedral") == 0) + mode = DIHEDRAL; + else if (strcmp(arg[3], "improper") == 0) + mode = IMPROPER; + else + error->all(FLERR, "Invalid compute count/type keyword {}", arg[3]); // error check - + if (mode == BOND && !atom->nbondtypes) - error->all(FLERR,"Compute count/type bond command with no bonds defined"); + error->all(FLERR, "Compute count/type bond command with no bonds defined"); if (mode == ANGLE && !atom->nangletypes) - error->all(FLERR,"Compute count/type bond command with no angles defined"); + error->all(FLERR, "Compute count/type bond command with no angles defined"); if (mode == DIHEDRAL && !atom->ndihedraltypes) - error->all(FLERR,"Compute count/type dihedral command with no dihedrals defined"); + error->all(FLERR, "Compute count/type dihedral command with no dihedrals defined"); if (mode == IMPROPER && !atom->nimpropertypes) - error->all(FLERR,"Compute count/type improper command with no impropers defined"); + error->all(FLERR, "Compute count/type improper command with no impropers defined"); // set vector lengths @@ -145,18 +151,23 @@ void ComputeCountType::compute_vector() int nvec; - if (mode == ATOM) nvec = count_atoms(); - else if (mode == BOND) nvec = count_bonds(); - else if (mode == ANGLE) nvec = count_angles(); - else if (mode == DIHEDRAL) nvec = count_dihedrals(); - else if (mode == IMPROPER) nvec = count_impropers(); - + if (mode == ATOM) + nvec = count_atoms(); + else if (mode == BOND) + nvec = count_bonds(); + else if (mode == ANGLE) + nvec = count_angles(); + else if (mode == DIHEDRAL) + nvec = count_dihedrals(); + else if (mode == IMPROPER) + nvec = count_impropers(); + // sum across procs as bigint, then convert to double // correct for multiple counting if newton_bond off for (int m = 0; m < nvec; m++) bcount_me[m] = count[m]; MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world); - + if (force->newton_bond == 0) { if (mode == BOND) for (int m = 0; m < nvec; m++) bcount[m] /= 2; @@ -165,7 +176,7 @@ void ComputeCountType::compute_vector() if (mode == DIHEDRAL || mode == IMPROPER) for (int m = 0; m < nvec; m++) bcount[m] /= 4; } - + for (int m = 0; m < nvec; m++) if (bcount[m] > MAXDOUBLEINT) error->all(FLERR, "Compute count/type overflow"); for (int m = 0; m < nvec; m++) vector[m] = bcount[m]; @@ -182,11 +193,10 @@ int ComputeCountType::count_atoms() int *mask = atom->mask; int nlocal = atom->nlocal; int ntypes = atom->ntypes; - + for (int m = 0; m < ntypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - count[type[i]-1]++; + if (mask[i] & groupbit) count[type[i] - 1]++; return ntypes; } @@ -206,34 +216,36 @@ int ComputeCountType::count_bonds() int *mask = atom->mask; int nlocal = atom->nlocal; int nbondtypes = atom->nbondtypes; - - int j,m,nbond,itype; + + int j, m, nbond, itype; int flag = 0; for (int m = 0; m < nbondtypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; for (m = 0; m < nbond; m++) { itype = bond_type[i][m]; if (itype == 0) continue; - + j = atom->map(bond_atom[i][m]); if (j < 0) { flag = 1; continue; } - + if ((mask[i] & groupbit) && (mask[j] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); - + if (flagany) error->all(FLERR, "Missing bond atom in compute count/type"); + return nbondtypes; } @@ -253,16 +265,16 @@ int ComputeCountType::count_angles() int *mask = atom->mask; int nlocal = atom->nlocal; int nangletypes = atom->nangletypes; - - int j1,j2,j3,m,nangle,itype; + + int j1, j2, j3, m, nangle, itype; int flag = 0; for (int m = 0; m < nangletypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { nangle = num_angle[i]; for (m = 0; m < nangle; m++) { itype = angle_type[i][m]; - + j1 = atom->map(angle_atom1[i][m]); j2 = atom->map(angle_atom2[i][m]); j3 = atom->map(angle_atom3[i][m]); @@ -270,19 +282,20 @@ int ComputeCountType::count_angles() flag = 1; continue; } - - if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && - (mask[j3] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing angle atom in compute count/type"); - + if (flagany) error->all(FLERR, "Missing angle atom in compute count/type"); + return nangletypes; } @@ -303,16 +316,16 @@ int ComputeCountType::count_dihedrals() int *mask = atom->mask; int nlocal = atom->nlocal; int ndihedraltypes = atom->ndihedraltypes; - - int j1,j2,j3,j4,m,ndihedral,itype; + + int j1, j2, j3, j4, m, ndihedral, itype; int flag = 0; for (int m = 0; m < ndihedraltypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { ndihedral = num_dihedral[i]; for (m = 0; m < ndihedral; m++) { itype = dihedral_type[i][m]; - + j1 = atom->map(dihedral_atom1[i][m]); j2 = atom->map(dihedral_atom2[i][m]); j3 = atom->map(dihedral_atom3[i][m]); @@ -321,19 +334,21 @@ int ComputeCountType::count_dihedrals() flag = 1; continue; } - - if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && - (mask[j3] & groupbit) && (mask[j4] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit) && + (mask[j4] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing dihedral atom in compute count/type"); - + if (flagany) error->all(FLERR, "Missing dihedral atom in compute count/type"); + return ndihedraltypes; } @@ -354,16 +369,16 @@ int ComputeCountType::count_impropers() int *mask = atom->mask; int nlocal = atom->nlocal; int nimpropertypes = atom->nimpropertypes; - - int j1,j2,j3,j4,m,nimproper,itype; + + int j1, j2, j3, j4, m, nimproper, itype; int flag = 0; for (int m = 0; m < nimpropertypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { nimproper = num_improper[i]; for (m = 0; m < nimproper; m++) { itype = improper_type[i][m]; - + j1 = atom->map(improper_atom1[i][m]); j2 = atom->map(improper_atom2[i][m]); j3 = atom->map(improper_atom3[i][m]); @@ -372,18 +387,20 @@ int ComputeCountType::count_impropers() flag = 1; continue; } - - if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && - (mask[j3] & groupbit) && (mask[j4] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit) && + (mask[j4] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing improper atom in compute count/type"); - + if (flagany) error->all(FLERR, "Missing improper atom in compute count/type"); + return nimpropertypes; } From e8a77c61ac72a9933f9d10d68026b5e0e9901214 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 16:19:53 -0400 Subject: [PATCH 33/41] whitespace --- doc/src/compute_count_type.rst | 2 +- src/lmptype.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index 2e3098b5ea..bb276ea21b 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -50,7 +50,7 @@ type: * :doc:`fix shake ` * :doc:`delete_bonds ` - + These commands can create and/or break topological bonds (angles, etc). In the case of breaking, they remove the bond (angle, etc) from the system, so that they no longer exist (:doc:`bond_style quartic diff --git a/src/lmptype.h b/src/lmptype.h index 891d5bdf89..ecd6ee5761 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -96,7 +96,7 @@ typedef int64_t bigint; #define MAXTAGINT INT_MAX #define MAXBIGINT INT64_MAX #define MAXDOUBLEINT 9007199254740992 // 2^53 - + #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT #define MPI_LMP_BIGINT MPI_LL From b6e211dd62147c4f35ee70396f2831dfc2613492 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 16:20:20 -0400 Subject: [PATCH 34/41] consistently declare variables when used only --- src/compute_count_type.cpp | 64 +++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 35 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 5c27dd39ca..47bbe36dd2 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -27,7 +27,7 @@ enum { ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER }; /* ---------------------------------------------------------------------- */ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), count(nullptr), bcount(nullptr), bcount_me(nullptr) + Compute(lmp, narg, arg), count(nullptr), bcount_me(nullptr), bcount(nullptr) { if (narg != 4) error->all(FLERR, "Incorrect number of args for compute count/type command"); @@ -117,24 +117,22 @@ double ComputeCountType::compute_scalar() int **bond_type = atom->bond_type; int nlocal = atom->nlocal; - int m, nbond; - int count = 0; - // count broken bonds with bond_type = 0 // ignore group setting since 2 atoms in a broken bond // can be arbitrarily far apart + int count = 0; for (int i = 0; i < nlocal; i++) { - nbond = num_bond[i]; - for (m = 0; m < nbond; m++) + int nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) if (bond_type[i][m] == 0) count++; } // sum across procs as bigint, then convert to double // correct for double counting if newton_bond off + bigint bcount = 0; bigint bcount_me = count; - bigint bcount; MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); if (force->newton_bond == 0) bcount /= 2; @@ -217,17 +215,16 @@ int ComputeCountType::count_bonds() int nlocal = atom->nlocal; int nbondtypes = atom->nbondtypes; - int j, m, nbond, itype; int flag = 0; for (int m = 0; m < nbondtypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { - nbond = num_bond[i]; - for (m = 0; m < nbond; m++) { - itype = bond_type[i][m]; + int nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) { + int itype = bond_type[i][m]; if (itype == 0) continue; - j = atom->map(bond_atom[i][m]); + int j = atom->map(bond_atom[i][m]); if (j < 0) { flag = 1; continue; @@ -266,18 +263,17 @@ int ComputeCountType::count_angles() int nlocal = atom->nlocal; int nangletypes = atom->nangletypes; - int j1, j2, j3, m, nangle, itype; int flag = 0; for (int m = 0; m < nangletypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { - nangle = num_angle[i]; - for (m = 0; m < nangle; m++) { - itype = angle_type[i][m]; + int nangle = num_angle[i]; + for (int m = 0; m < nangle; m++) { + int itype = angle_type[i][m]; - j1 = atom->map(angle_atom1[i][m]); - j2 = atom->map(angle_atom2[i][m]); - j3 = atom->map(angle_atom3[i][m]); + int j1 = atom->map(angle_atom1[i][m]); + int j2 = atom->map(angle_atom2[i][m]); + int j3 = atom->map(angle_atom3[i][m]); if (j1 < 0 || j2 < 0 || j3 < 0) { flag = 1; continue; @@ -317,19 +313,18 @@ int ComputeCountType::count_dihedrals() int nlocal = atom->nlocal; int ndihedraltypes = atom->ndihedraltypes; - int j1, j2, j3, j4, m, ndihedral, itype; int flag = 0; for (int m = 0; m < ndihedraltypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { - ndihedral = num_dihedral[i]; - for (m = 0; m < ndihedral; m++) { - itype = dihedral_type[i][m]; + int ndihedral = num_dihedral[i]; + for (int m = 0; m < ndihedral; m++) { + int itype = dihedral_type[i][m]; - j1 = atom->map(dihedral_atom1[i][m]); - j2 = atom->map(dihedral_atom2[i][m]); - j3 = atom->map(dihedral_atom3[i][m]); - j4 = atom->map(dihedral_atom4[i][m]); + int j1 = atom->map(dihedral_atom1[i][m]); + int j2 = atom->map(dihedral_atom2[i][m]); + int j3 = atom->map(dihedral_atom3[i][m]); + int j4 = atom->map(dihedral_atom4[i][m]); if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) { flag = 1; continue; @@ -370,19 +365,18 @@ int ComputeCountType::count_impropers() int nlocal = atom->nlocal; int nimpropertypes = atom->nimpropertypes; - int j1, j2, j3, j4, m, nimproper, itype; int flag = 0; for (int m = 0; m < nimpropertypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { - nimproper = num_improper[i]; - for (m = 0; m < nimproper; m++) { - itype = improper_type[i][m]; + int nimproper = num_improper[i]; + for (int m = 0; m < nimproper; m++) { + int itype = improper_type[i][m]; - j1 = atom->map(improper_atom1[i][m]); - j2 = atom->map(improper_atom2[i][m]); - j3 = atom->map(improper_atom3[i][m]); - j4 = atom->map(improper_atom4[i][m]); + int j1 = atom->map(improper_atom1[i][m]); + int j2 = atom->map(improper_atom2[i][m]); + int j3 = atom->map(improper_atom3[i][m]); + int j4 = atom->map(improper_atom4[i][m]); if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) { flag = 1; continue; From 3f8cd4577cd161a7ac9a04781bcf2882336a46fb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 16:20:29 -0400 Subject: [PATCH 35/41] update unit test --- unittest/commands/test_compute_global.cpp | 66 ++++++++++++++++++----- 1 file changed, 52 insertions(+), 14 deletions(-) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index 7e5b3d4483..6c365c8c2b 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -308,45 +308,83 @@ TEST_F(ComputeGlobalTest, Counts) command("variable t3 atom type==3"); command("variable t4 atom type==4"); command("variable t5 atom type==5"); - command("compute asum all reduce sum v_t1 v_t2 v_t3 v_t4 v_t5"); - command("compute acnt all count/type atom"); + command("compute tsum all reduce sum v_t1 v_t2 v_t3 v_t4 v_t5"); + command("compute tcnt all count/type atom"); command("compute bcnt all count/type bond"); - command("thermo_style custom c_asum[*] c_acnt[*]"); + command("compute acnt all count/type angle"); + command("compute dcnt all count/type dihedral"); + command("compute icnt all count/type improper"); + command("thermo_style custom c_tsum[*] c_tcnt[*] c_bcnt[*] c_acnt[*] c_dcnt[*] c_icnt[*]"); command("run 0 post no"); END_HIDE_OUTPUT(); - auto asum = get_vector("asum"); - auto acnt = get_vector("acnt"); + auto tsum = get_vector("tsum"); + auto tcnt = get_vector("tcnt"); auto bcnt = get_vector("bcnt"); auto bbrk = get_scalar("bcnt"); + auto acnt = get_vector("acnt"); + auto dcnt = get_vector("dcnt"); + auto icnt = get_vector("icnt"); + + EXPECT_DOUBLE_EQ(tsum[0], tcnt[0]); + EXPECT_DOUBLE_EQ(tsum[1], tcnt[1]); + EXPECT_DOUBLE_EQ(tsum[2], tcnt[2]); + EXPECT_DOUBLE_EQ(tsum[3], tcnt[3]); + EXPECT_DOUBLE_EQ(tsum[4], tcnt[4]); EXPECT_DOUBLE_EQ(bbrk, 0.0); - EXPECT_DOUBLE_EQ(asum[0], acnt[0]); - EXPECT_DOUBLE_EQ(asum[1], acnt[1]); - EXPECT_DOUBLE_EQ(asum[2], acnt[2]); - EXPECT_DOUBLE_EQ(asum[3], acnt[3]); - EXPECT_DOUBLE_EQ(asum[4], acnt[4]); - EXPECT_DOUBLE_EQ(bcnt[0], 3.0); EXPECT_DOUBLE_EQ(bcnt[1], 6.0); EXPECT_DOUBLE_EQ(bcnt[2], 3.0); EXPECT_DOUBLE_EQ(bcnt[3], 2.0); EXPECT_DOUBLE_EQ(bcnt[4], 10.0); + EXPECT_DOUBLE_EQ(acnt[0], 6.0); + EXPECT_DOUBLE_EQ(acnt[1], 10.0); + EXPECT_DOUBLE_EQ(acnt[2], 5.0); + EXPECT_DOUBLE_EQ(acnt[3], 9.0); + + EXPECT_DOUBLE_EQ(dcnt[0], 3.0); + EXPECT_DOUBLE_EQ(dcnt[1], 8.0); + EXPECT_DOUBLE_EQ(dcnt[2], 3.0); + EXPECT_DOUBLE_EQ(dcnt[3], 4.0); + EXPECT_DOUBLE_EQ(dcnt[4], 13.0); + + EXPECT_DOUBLE_EQ(icnt[0], 1.0); + EXPECT_DOUBLE_EQ(icnt[1], 1.0); + BEGIN_HIDE_OUTPUT(); - command("delete_bonds all bond 3"); + command("delete_bonds all bond 3 remove"); command("run 0 post no"); END_HIDE_OUTPUT(); bcnt = get_vector("bcnt"); bbrk = get_scalar("bcnt"); - EXPECT_DOUBLE_EQ(bbrk, 0.0); // should be 3 + acnt = get_vector("acnt"); + dcnt = get_vector("dcnt"); + icnt = get_vector("icnt"); + + EXPECT_DOUBLE_EQ(bbrk, 0.0); EXPECT_DOUBLE_EQ(bcnt[0], 3.0); EXPECT_DOUBLE_EQ(bcnt[1], 6.0); - EXPECT_DOUBLE_EQ(bcnt[2], 3.0); // should be 0 + EXPECT_DOUBLE_EQ(bcnt[2], 0.0); EXPECT_DOUBLE_EQ(bcnt[3], 2.0); EXPECT_DOUBLE_EQ(bcnt[4], 10.0); + + EXPECT_DOUBLE_EQ(acnt[0], 6.0); + EXPECT_DOUBLE_EQ(acnt[1], 10.0); + EXPECT_DOUBLE_EQ(acnt[2], 5.0); + EXPECT_DOUBLE_EQ(acnt[3], 9.0); + + EXPECT_DOUBLE_EQ(dcnt[0], 3.0); + EXPECT_DOUBLE_EQ(dcnt[1], 8.0); + EXPECT_DOUBLE_EQ(dcnt[2], 3.0); + EXPECT_DOUBLE_EQ(dcnt[3], 4.0); + EXPECT_DOUBLE_EQ(dcnt[4], 13.0); + + EXPECT_DOUBLE_EQ(icnt[0], 1.0); + EXPECT_DOUBLE_EQ(icnt[1], 1.0); } } // namespace LAMMPS_NS From 6d050374c34b18ea65b06fe06feed36834d61e35 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 16:24:43 -0400 Subject: [PATCH 36/41] fix compilation with -DLAMMPS_BIGBIG --- src/compute_count_type.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 47bbe36dd2..84b31e8b85 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -208,7 +208,7 @@ int ComputeCountType::count_atoms() int ComputeCountType::count_bonds() { - int **bond_atom = atom->bond_atom; + tagint **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; int *num_bond = atom->num_bond; int *mask = atom->mask; @@ -254,9 +254,9 @@ int ComputeCountType::count_bonds() int ComputeCountType::count_angles() { - int **angle_atom1 = atom->angle_atom1; - int **angle_atom2 = atom->angle_atom2; - int **angle_atom3 = atom->angle_atom3; + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom2 = atom->angle_atom2; + tagint **angle_atom3 = atom->angle_atom3; int **angle_type = atom->angle_type; int *num_angle = atom->num_angle; int *mask = atom->mask; @@ -303,10 +303,10 @@ int ComputeCountType::count_angles() int ComputeCountType::count_dihedrals() { - int **dihedral_atom1 = atom->dihedral_atom1; - int **dihedral_atom2 = atom->dihedral_atom2; - int **dihedral_atom3 = atom->dihedral_atom3; - int **dihedral_atom4 = atom->dihedral_atom4; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; int **dihedral_type = atom->dihedral_type; int *num_dihedral = atom->num_dihedral; int *mask = atom->mask; @@ -355,10 +355,10 @@ int ComputeCountType::count_dihedrals() int ComputeCountType::count_impropers() { - int **improper_atom1 = atom->improper_atom1; - int **improper_atom2 = atom->improper_atom2; - int **improper_atom3 = atom->improper_atom3; - int **improper_atom4 = atom->improper_atom4; + tagint **improper_atom1 = atom->improper_atom1; + tagint **improper_atom2 = atom->improper_atom2; + tagint **improper_atom3 = atom->improper_atom3; + tagint **improper_atom4 = atom->improper_atom4; int **improper_type = atom->improper_type; int *num_improper = atom->num_improper; int *mask = atom->mask; From 36632e33566bd6a93f08539461d5051bd71a1097 Mon Sep 17 00:00:00 2001 From: Joel Thomas Clemmer Date: Tue, 2 May 2023 14:28:17 -0600 Subject: [PATCH 37/41] Clarifying how BPM package works --- doc/src/Howto_broken_bonds.rst | 14 ++++++++------ doc/src/compute_count_type.rst | 14 ++++++++------ 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/doc/src/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index 282d2b32a4..9fcd676c36 100644 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -19,6 +19,14 @@ in the broken bond are "turned on". Angles, dihedrals, etc cannot be defined for a system when :doc:`bond_style quartic ` is used. +Similarly, bond styles in the BPM package are also incompatible with +angles, dihedrals, etc. and when a bond breaks its type is set to zero. +However, in the BPM package one can either turn off all pair interactions +between bonded particles or leave them on, overlaying pair forces on +top of bond forces. To remove pair forces, the special bond list is +dynamically updated. More details can be found on the :doc:`Howto BPM +` page. + The :doc:`fix bond/break ` and :doc:`fix bond/react ` commands allow breaking of bonds within a molecular topology with may also define angles, dihedrals, etc. These commands @@ -27,12 +35,6 @@ well as the appropriate angle, dihederal, etc interactions which include the bond. They also trigger a rebuild of the neighbor list when this occurs, to turn on the appropriate pairwise forces. -In the BPM package, one can either turn off all pair interactions -between bonded particles or leave them on, overlaying pair forces on -top of bond forces. To remove pair forces, the special bond list is -dynamically updated. More details can be found on the :doc:`Howto BPM -` page. - Note that when bonds are dumped to a file via the :doc:`dump local ` command, bonds with type 0 are not included. diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index bb276ea21b..b92fe2f2fb 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -54,8 +54,9 @@ type: These commands can create and/or break topological bonds (angles, etc). In the case of breaking, they remove the bond (angle, etc) from the system, so that they no longer exist (:doc:`bond_style quartic -` is the exception, see the discussion below). Thus -they are not included in the counts for each type: +` and :doc:`BPM bond styles ` are exceptions, +see the discussion below). Thus they are not included in the counts +for each type: * :doc:`delete_bonds remove ` * :doc:`bond_style quartic ` @@ -74,10 +75,11 @@ type is tallied. Only bonds with both atoms in the specified group are counted. For {mode} = {bond}, broken bonds with a bond type of zero are also -counted. The :doc:`bond_style quartic ` breaks a bond -by doing this. See the :doc:`Howto broken bonds ` -doc page for more details. Note that the group setting is ignored for -broken bonds; all broken bonds in the system are counted. +counted. The :doc:`bond_style quartic ` and :doc:`BPM +bond styles ` break bonds by doing this. See the :doc:` +Howto broken bonds ` doc page for more details. +Note that the group setting is ignored for broken bonds; all broken +bonds in the system are counted. If the {mode} setting is {angle} then the count of angles for each angle type is tallied. Only angles with all 3 atoms in the specified From 7b5fecfad72d6313536edc4d6c3a0a8595f37b15 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 17:24:38 -0400 Subject: [PATCH 38/41] make phana installation in LAMMPS compatible with ubuntu 18.04LTS --- cmake/Modules/Tools.cmake | 4 ++-- tools/phonon/CMakeLists.spglib | 2 +- tools/phonon/CMakeLists.txt | 2 +- tools/phonon/tricubic/CMakeLists.txt | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/cmake/Modules/Tools.cmake b/cmake/Modules/Tools.cmake index ba91a557bc..285c5f2405 100644 --- a/cmake/Modules/Tools.cmake +++ b/cmake/Modules/Tools.cmake @@ -33,6 +33,8 @@ if(BUILD_TOOLS) endif() install(TARGETS msi2lmp DESTINATION ${CMAKE_INSTALL_BINDIR}) install(FILES ${LAMMPS_DOC_DIR}/msi2lmp.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1) + + add_subdirectory(${LAMMPS_TOOLS_DIR}/phonon ${CMAKE_BINARY_DIR}/phana_build) endif() if(BUILD_LAMMPS_SHELL) @@ -63,8 +65,6 @@ if(BUILD_LAMMPS_SHELL) install(TARGETS lammps-shell EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR}) install(DIRECTORY ${LAMMPS_TOOLS_DIR}/lammps-shell/icons DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/) install(FILES ${LAMMPS_TOOLS_DIR}/lammps-shell/lammps-shell.desktop DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/applications/) - - add_subdirectory(${LAMMPS_TOOLS_DIR}/phonon ${CMAKE_BINARY_DIR}/phana_build) endif() diff --git a/tools/phonon/CMakeLists.spglib b/tools/phonon/CMakeLists.spglib index 9e22445b80..566e58a2b6 100644 --- a/tools/phonon/CMakeLists.spglib +++ b/tools/phonon/CMakeLists.spglib @@ -47,7 +47,7 @@ set(SOURCES ${PROJECT_SOURCE_DIR}/src/arithmetic.c ${PROJECT_SOURCE_DIR}/src/symmetry.c) add_library(symspg STATIC ${SOURCES}) -install(TARGETS symspg LIBRARY) +install(TARGETS symspg ARCHIVE DESTINATION ${CMAKE_INSTALL_PREFIX}/lib) # Header file install(FILES ${PROJECT_SOURCE_DIR}/src/spglib.h DESTINATION ${CMAKE_INSTALL_PREFIX}/include) diff --git a/tools/phonon/CMakeLists.txt b/tools/phonon/CMakeLists.txt index 02ccee25a3..60da1cc79f 100644 --- a/tools/phonon/CMakeLists.txt +++ b/tools/phonon/CMakeLists.txt @@ -1,7 +1,7 @@ # Support Linux from Ubuntu 20.04LTS onward, CentOS 7.x (with EPEL), # macOS, MSVC 2019 (=Version 16) -cmake_minimum_required(VERSION 3.16) +cmake_minimum_required(VERSION 3.10) # set timestamp of downloaded files to that of archive if(POLICY CMP0135) diff --git a/tools/phonon/tricubic/CMakeLists.txt b/tools/phonon/tricubic/CMakeLists.txt index d70b44c2b2..190c6ceb80 100644 --- a/tools/phonon/tricubic/CMakeLists.txt +++ b/tools/phonon/tricubic/CMakeLists.txt @@ -1,7 +1,7 @@ # Support Linux from Ubuntu 20.04LTS onward, CentOS 7.x (with EPEL), # macOS, MSVC 2019 (=Version 16) -cmake_minimum_required(VERSION 3.16) +cmake_minimum_required(VERSION 3.10) # set up project project(tricubic VERSION 1.1 DESCRIPTION "Tricubic library" LANGUAGES CXX) From a8dcb3987a2279db88cfa3d684522ef1912b5ff1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 17:29:06 -0400 Subject: [PATCH 39/41] silence compiler warning --- tools/phonon/phonopy.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tools/phonon/phonopy.cpp b/tools/phonon/phonopy.cpp index a6b991efac..2cee319aa7 100644 --- a/tools/phonon/phonopy.cpp +++ b/tools/phonon/phonopy.cpp @@ -250,7 +250,8 @@ void Phonopy::phonopy() for (int jdim = 0; jdim < 3; ++jdim) fprintf(fp, "%lg ", dm->basevec[ndim++]); fprintf(fp, "\n"); } - for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]); fprintf(fp, "\n"); + for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]); + fprintf(fp, "\n"); for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "%d ", num_type[ip]); fprintf(fp, "\nDirect\n"); for (int ip = 0; ip < ntype; ++ip){ From 3e3ed89f33428e6324bca28f90b3bc966f7c804a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 17:30:23 -0400 Subject: [PATCH 40/41] avoid out-of-range access if an angle/dihedral/improper type is set to 0 --- src/compute_count_type.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 84b31e8b85..3cb2e6bc3a 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -282,7 +282,7 @@ int ComputeCountType::count_angles() if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit)) { if (itype > 0) count[itype - 1]++; - else + else if (itype < 0) count[-itype - 1]++; } } @@ -334,7 +334,7 @@ int ComputeCountType::count_dihedrals() (mask[j4] & groupbit)) { if (itype > 0) count[itype - 1]++; - else + else if (itype < 0) count[-itype - 1]++; } } @@ -386,7 +386,7 @@ int ComputeCountType::count_impropers() (mask[j4] & groupbit)) { if (itype > 0) count[itype - 1]++; - else + else if (itype < 0) count[-itype - 1]++; } } From 7ab30aa4684bf2bbf843f1b930274f5d0aa500de Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 17:32:57 -0400 Subject: [PATCH 41/41] add versionadded tag --- doc/src/compute_count_type.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index b92fe2f2fb..73e9d50f73 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -25,6 +25,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + Define a computation that counts the current number of atoms for each atom type. Or the number of bonds (angles, dihedrals, impropers) for each bond (angle, dihedral, improper) type.