From cae83c25d9bf9fe59e45262c20339b35748687b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kliment=20Olechnovi=C4=8D?= <54395582+kliment-olechnovic@users.noreply.github.com> Date: Mon, 7 Oct 2024 11:14:03 +0200 Subject: [PATCH] Update Voronota-LT, add alternative modes to Voronoi analysis. (#456) Version 0.9.4 of Voronota-LT has an improved API and a new ability to update the tessellation without fully rebuilding it when possible. The Voronoi analysis class now has adjustable running mode. The class implementation makes use of the "pImpl" idiom. Examples collection now contains "voronota" directory with polymer-based and peptide based multichain simulation configuration files and the corresponding output logs. --- docs/_docs/analysis.md | 1 + docs/schema.yml | 5 + examples/voronota/voronota_peptides.out.json | 2696 +++++++++++++++++ examples/voronota/voronota_peptides.yml | 57 + examples/voronota/voronota_polymers.out.json | 1396 +++++++++ examples/voronota/voronota_polymers.yml | 83 + src/analysis.h | 29 +- src/voronota.cpp | 265 +- src/voronotalt/VERSION.txt | 2 +- src/voronotalt/basic_types_and_functions.h | 112 +- src/voronotalt/conversion_to_input.h | 12 +- .../parallelization_configuration.h | 8 + src/voronotalt/periodic_box.h | 94 + src/voronotalt/preparation_for_tessellation.h | 193 -- src/voronotalt/radical_tessellation.h | 348 ++- ...adical_tessellation_contact_construction.h | 69 +- src/voronotalt/simplified_aw_tessellation.h | 48 +- ...ied_aw_tessellation_contact_construction.h | 58 +- src/voronotalt/spheres_container.h | 664 ++++ src/voronotalt/spheres_searcher.h | 293 +- src/voronotalt/time_recorder.h | 6 +- .../updateable_radical_tessellation.h | 671 ++++ src/voronotalt/voronotalt.h | 1 + 23 files changed, 6583 insertions(+), 528 deletions(-) create mode 100644 examples/voronota/voronota_peptides.out.json create mode 100644 examples/voronota/voronota_peptides.yml create mode 100644 examples/voronota/voronota_polymers.out.json create mode 100755 examples/voronota/voronota_polymers.yml create mode 100644 src/voronotalt/parallelization_configuration.h create mode 100644 src/voronotalt/periodic_box.h delete mode 100644 src/voronotalt/preparation_for_tessellation.h create mode 100644 src/voronotalt/spheres_container.h create mode 100644 src/voronotalt/updateable_radical_tessellation.h diff --git a/docs/_docs/analysis.md b/docs/_docs/analysis.md index 362032238..80a956d74 100644 --- a/docs/_docs/analysis.md +++ b/docs/_docs/analysis.md @@ -377,6 +377,7 @@ The algorithm is significantly faster than the above `sasa` analysis. `nstep` | Interval between samples `nskip=0` | Number of initial steps excluded from the analysis `file` | Optionally stream surface area for each `nstep` to file (`.dat|.dat.gz`) +`mode=full` | Running mode: `full`, `interchain`, `updateable` `radius=1.4` | Probe radius (Å) diff --git a/docs/schema.yml b/docs/schema.yml index 500f5d65b..66084c266 100644 --- a/docs/schema.yml +++ b/docs/schema.yml @@ -1312,6 +1312,11 @@ properties: type: string pattern: "(.*?)\\.(dat|dat.gz)$" description: "Streaming filename (.dat|.dat.gz)" + mode: + type: string + enum: ["full", "interchain", "updateable"] + default: "full" + description: "Running mode" required: [nstep] additionalProperties: false diff --git a/examples/voronota/voronota_peptides.out.json b/examples/voronota/voronota_peptides.out.json new file mode 100644 index 000000000..e531e8fbf --- /dev/null +++ b/examples/voronota/voronota_peptides.out.json @@ -0,0 +1,2696 @@ +{ + "analysis": [ + { + "sanity": { + "nstep": 10, + "samples": 1000 + } + }, + { + "xtcfile": { + "file": "peptides_traj.xtc", + "nstep": 1, + "samples": 10000 + } + }, + { + "savestate": { + "file": "peptides_structure.gro" + } + }, + { + "sasa": { + "molecule": "peptide", + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/dbjh", + "relative time": 0.0645, + "samples": 10000, + "slices_per_atom": 20, + "⟨SASA²⟩-⟨SASA⟩²": 123516.51044967347, + "⟨SASA⟩": 10739.539171115213 + } + }, + { + "voronota": { + "mode": "full", + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/mq8k", + "relative time": 0.0206, + "samples": 10000, + "⟨area_of_contacts⟩": 10482.398435376826, + "⟨area_of_sas_per_chain⟩": 10739.479134995941, + "⟨area_of_sas⟩": 53697.39567497984, + "⟨volume_inside_sas⟩": 114871.95801628771 + } + }, + { + "voronota": { + "mode": "interchain", + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/mq8k", + "relative time": 0.0135, + "samples": 10000, + "⟨area_of_contacts⟩": 4.097313868222107, + "⟨area_of_sas_per_chain⟩": null, + "⟨area_of_sas⟩": null, + "⟨volume_inside_sas⟩": null + } + }, + { + "voronota": { + "mode": "updateable", + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/mq8k", + "relative time": 0.0212, + "samples": 10000, + "⟨area_of_contacts⟩": 10482.398435376826, + "⟨area_of_sas_per_chain⟩": 10739.479134995941, + "⟨area_of_sas⟩": 53697.39567497984, + "⟨volume_inside_sas⟩": 114871.95801628771 + } + } + ], + "compiler": "13.2.0", + "energy": [ + { + "hamiltonian": [ + { + "bonded": { + "bondlist-intramolecular": [ + { + "harmonic": { + "index": [ + 0, + 1 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 1, + 2 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 2, + 3 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 3, + 4 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 4, + 5 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 5, + 6 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 6, + 7 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 7, + 8 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 8, + 9 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 9, + 10 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 10, + 11 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 11, + 12 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 12, + 13 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 13, + 14 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 14, + 15 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 15, + 16 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 16, + 17 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 17, + 18 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 18, + 19 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 19, + 20 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 20, + 21 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 21, + 22 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 22, + 23 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 23, + 24 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 24, + 25 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 25, + 26 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 26, + 27 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 27, + 28 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 28, + 29 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 29, + 30 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 30, + 31 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 31, + 32 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 32, + 33 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 33, + 34 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 34, + 35 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 35, + 36 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 36, + 37 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 37, + 38 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 38, + 39 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 39, + 40 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 40, + 41 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 41, + 42 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 42, + 43 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 43, + 44 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 44, + 45 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 45, + 46 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 46, + 47 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 47, + 48 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 48, + 49 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 50, + 51 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 51, + 52 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 52, + 53 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 53, + 54 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 54, + 55 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 55, + 56 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 56, + 57 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 57, + 58 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 58, + 59 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 59, + 60 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 60, + 61 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 61, + 62 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 62, + 63 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 63, + 64 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 64, + 65 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 65, + 66 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 66, + 67 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 67, + 68 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 68, + 69 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 69, + 70 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 70, + 71 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 71, + 72 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 72, + 73 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 73, + 74 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 74, + 75 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 75, + 76 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 76, + 77 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 77, + 78 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 78, + 79 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 79, + 80 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 80, + 81 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 81, + 82 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 82, + 83 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 83, + 84 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 84, + 85 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 85, + 86 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 86, + 87 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 87, + 88 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 88, + 89 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 89, + 90 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 90, + 91 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 91, + 92 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 92, + 93 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 93, + 94 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 94, + 95 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 95, + 96 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 96, + 97 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 97, + 98 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 98, + 99 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 100, + 101 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 101, + 102 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 102, + 103 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 103, + 104 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 104, + 105 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 105, + 106 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 106, + 107 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 107, + 108 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 108, + 109 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 109, + 110 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 110, + 111 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 111, + 112 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 112, + 113 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 113, + 114 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 114, + 115 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 115, + 116 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 116, + 117 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 117, + 118 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 118, + 119 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 119, + 120 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 120, + 121 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 121, + 122 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 122, + 123 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 123, + 124 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 124, + 125 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 125, + 126 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 126, + 127 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 127, + 128 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 128, + 129 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 129, + 130 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 130, + 131 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 131, + 132 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 132, + 133 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 133, + 134 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 134, + 135 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 135, + 136 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 136, + 137 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 137, + 138 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 138, + 139 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 139, + 140 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 140, + 141 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 141, + 142 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 142, + 143 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 143, + 144 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 144, + 145 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 145, + 146 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 146, + 147 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 147, + 148 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 148, + 149 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 150, + 151 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 151, + 152 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 152, + 153 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 153, + 154 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 154, + 155 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 155, + 156 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 156, + 157 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 157, + 158 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 158, + 159 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 159, + 160 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 160, + 161 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 161, + 162 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 162, + 163 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 163, + 164 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 164, + 165 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 165, + 166 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 166, + 167 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 167, + 168 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 168, + 169 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 169, + 170 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 170, + 171 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 171, + 172 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 172, + 173 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 173, + 174 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 174, + 175 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 175, + 176 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 176, + 177 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 177, + 178 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 178, + 179 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 179, + 180 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 180, + 181 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 181, + 182 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 182, + 183 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 183, + 184 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 184, + 185 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 185, + 186 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 186, + 187 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 187, + 188 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 188, + 189 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 189, + 190 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 190, + 191 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 191, + 192 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 192, + 193 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 193, + 194 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 194, + 195 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 195, + 196 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 196, + 197 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 197, + 198 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 198, + 199 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 200, + 201 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 201, + 202 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 202, + 203 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 203, + 204 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 204, + 205 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 205, + 206 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 206, + 207 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 207, + 208 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 208, + 209 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 209, + 210 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 210, + 211 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 211, + 212 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 212, + 213 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 213, + 214 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 214, + 215 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 215, + 216 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 216, + 217 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 217, + 218 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 218, + 219 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 219, + 220 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 220, + 221 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 221, + 222 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 222, + 223 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 223, + 224 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 224, + 225 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 225, + 226 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 226, + 227 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 227, + 228 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 228, + 229 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 229, + 230 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 230, + 231 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 231, + 232 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 232, + 233 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 233, + 234 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 234, + 235 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 235, + 236 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 236, + 237 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 237, + 238 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 238, + 239 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 239, + 240 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 240, + 241 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 241, + 242 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 242, + 243 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 243, + 244 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 244, + 245 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 245, + 246 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 246, + 247 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 247, + 248 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + }, + { + "harmonic": { + "index": [ + 248, + 249 + ], + "k": 3.0000000000000004, + "req": 7.0 + } + } + ], + "relative time": 0.002811205095690397 + } + }, + { + "nonbonded": { + "functor potential": { + "default": [ + { + "wca": { + "mixing": "LB" + } + } + ] + }, + "relative time": 0.4171505696957756, + "selfenergy": { + "dipole": false, + "monopole": false + }, + "summation_policy": "serial" + } + } + ] + } + ], + "geometry": { + "length": [ + 300.0, + 350.0, + 400.0 + ], + "type": "cuboid" + }, + "git revision": "3d460516 (2024-07-24)", + "groups": [ + { + "peptide": { + "compressible": false, + "index": [ + 0, + 49 + ], + "size": 50 + } + }, + { + "peptide": { + "compressible": false, + "index": [ + 50, + 99 + ], + "size": 50 + } + }, + { + "peptide": { + "compressible": false, + "index": [ + 100, + 149 + ], + "size": 50 + } + }, + { + "peptide": { + "compressible": false, + "index": [ + 150, + 199 + ], + "size": 50 + } + }, + { + "peptide": { + "compressible": false, + "index": [ + 200, + 249 + ], + "size": 50 + } + } + ], + "montecarlo": { + "average potential energy (kT)": 17638.4008666005, + "last move": "transrot" + }, + "moves": [ + { + "transrot": { + "acceptance": 0.637, + "dir": [ + 1.0, + 1.0, + 1.0 + ], + "dp": 0.0, + "dprot": 0.0, + "molecule": "peptide", + "molid": 0, + "moves": 2500545, + "relative time": 0.0435, + "repeat": 250, + "stochastic": true, + "√⟨r²⟩": 0.795 + } + }, + { + "pivot": { + "acceptance": 0.953, + "dprot": 0.5, + "molecule": "peptide", + "moves": 2449471, + "relative time": 0.779, + "repeat": 245, + "skipped": 2.29e-05, + "stochastic": true, + "√⟨r_cm²⟩": 2.71 + } + }, + { + "moltransrot": { + "R/Å ≈ √(⟨r²⟩/4⟨θ²⟩)": 40.0, + "acceptance": 0.985, + "dir": [ + 1.0, + 1.0, + 1.0 + ], + "dirrot": [ + 0.0, + 0.0, + 0.0 + ], + "dp": 20.0, + "dprot": 0.5, + "molecule": "peptide", + "molid": 0, + "moves": 99984, + "relative time": 0.0523, + "repeat": 10, + "stochastic": true, + "√⟨r²⟩": 11.4, + "√⟨θ²⟩": 0.143, + "√⟨θ²⟩/°": 8.2 + } + } + ], + "number of groups": 5, + "number of particles": 250, + "number of sweeps": 10000, + "reactionlist": null, + "relative drift": -7.076029808911516e-14, + "simulation time": { + "in minutes": 12.3, + "in seconds": 738 + }, + "temperature": 298.0 +} diff --git a/examples/voronota/voronota_peptides.yml b/examples/voronota/voronota_peptides.yml new file mode 100644 index 000000000..952c6f3ce --- /dev/null +++ b/examples/voronota/voronota_peptides.yml @@ -0,0 +1,57 @@ +#!/usr/bin/env yason.py +{% set N = 50 %} + +temperature: 298 +random: {seed: hardware} +geometry: {type: cuboid, length: [300, 350, 400]} +mcloop: {macro: 10, micro: 1000} + +atomlist: + - GLN : { sigma: 6.0, mw: 54.0, dp: 2.0, eps: 1.0 } + +moleculelist: + - peptide: + structure: + fasta: {% for i in range(N) %}Q{% endfor %} + k: 3.0 + req: 7.0 + +insertmolecules: + - peptide: {N: 5} + +energy: + - bonded: {} + - nonbonded: + default: + - wca: {mixing: LB} + +moves: + - transrot: {molecule: peptide, repeat: N} + - pivot: {molecule: peptide, dprot: 0.5, repeat: N} + - moltransrot: {molecule: peptide, dp: 20.0, dprot: 0.5, repeat: 10} + +analysis: + - sanity: {nstep: 10} + - xtcfile: {file: peptides_traj.xtc, nstep: 1} + - savestate: {file: peptides_structure.gro} + - sasa: + file: peptides_sasa.dat + policy: molecular + molecule: peptide + radius: 2.0 + nstep: 1 + - voronoi: + file: peptides_voronoi_full.dat + mode: full + radius: 2.0 + nstep: 1 + - voronoi: + file: peptides_voronoi_interchain.dat + mode: interchain + radius: 2.0 + nstep: 1 + - voronoi: + file: peptides_voronoi_updateable.dat + mode: updateable + radius: 2.0 + nstep: 1 diff --git a/examples/voronota/voronota_polymers.out.json b/examples/voronota/voronota_polymers.out.json new file mode 100644 index 000000000..9020d6959 --- /dev/null +++ b/examples/voronota/voronota_polymers.out.json @@ -0,0 +1,1396 @@ +{ + "analysis": [ + { + "sanity": { + "nstep": 10, + "samples": 1000 + } + }, + { + "xtcfile": { + "file": "polymers_traj.xtc", + "nstep": 1, + "samples": 10000 + } + }, + { + "savestate": { + "file": "polymers_structure.pqr" + } + }, + { + "sasa": { + "molecule": "polymer", + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/dbjh", + "relative time": 0.3, + "samples": 10000, + "slices_per_atom": 20, + "⟨SASA²⟩-⟨SASA⟩²": 76035.48336485625, + "⟨SASA⟩": 3668.374115022107 + } + }, + { + "voronota": { + "mode": "full", + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/mq8k", + "relative time": 0.027, + "samples": 10000, + "⟨area_of_contacts⟩": 3502.777808040449, + "⟨area_of_sas_per_chain⟩": 3668.364917632628, + "⟨area_of_sas⟩": 29346.919341061024, + "⟨volume_inside_sas⟩": 57270.355681232155 + } + }, + { + "voronota": { + "mode": "interchain", + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/mq8k", + "relative time": 0.0102, + "samples": 10000, + "⟨area_of_contacts⟩": 138.2561356671105, + "⟨area_of_sas_per_chain⟩": null, + "⟨area_of_sas⟩": null, + "⟨volume_inside_sas⟩": null + } + }, + { + "voronota": { + "mode": "updateable", + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/mq8k", + "relative time": 0.0242, + "samples": 10000, + "⟨area_of_contacts⟩": 3502.777808040449, + "⟨area_of_sas_per_chain⟩": 3668.364917632628, + "⟨area_of_sas⟩": 29346.919341061024, + "⟨volume_inside_sas⟩": 57270.355681232155 + } + } + ], + "compiler": "13.2.0", + "energy": [ + { + "hamiltonian": [ + { + "ContainerOverlap": { + "relative time": 1.5625280024417666e-05 + } + }, + { + "bonded": { + "bondlist-intramolecular": [ + { + "harmonic": { + "index": [ + 0, + 1 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 1, + 2 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 2, + 3 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 3, + 4 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 4, + 5 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 5, + 6 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 6, + 7 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 7, + 8 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 8, + 9 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 9, + 10 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 10, + 11 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 11, + 12 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 12, + 13 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 13, + 14 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 15, + 16 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 16, + 17 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 17, + 18 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 18, + 19 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 19, + 20 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 20, + 21 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 21, + 22 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 22, + 23 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 23, + 24 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 24, + 25 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 25, + 26 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 26, + 27 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 27, + 28 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 28, + 29 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 30, + 31 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 31, + 32 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 32, + 33 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 33, + 34 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 34, + 35 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 35, + 36 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 36, + 37 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 37, + 38 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 38, + 39 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 39, + 40 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 40, + 41 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 41, + 42 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 42, + 43 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 43, + 44 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 45, + 46 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 46, + 47 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 47, + 48 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 48, + 49 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 49, + 50 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 50, + 51 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 51, + 52 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 52, + 53 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 53, + 54 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 54, + 55 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 55, + 56 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 56, + 57 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 57, + 58 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 58, + 59 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 60, + 61 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 61, + 62 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 62, + 63 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 63, + 64 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 64, + 65 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 65, + 66 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 66, + 67 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 67, + 68 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 68, + 69 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 69, + 70 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 70, + 71 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 71, + 72 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 72, + 73 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 73, + 74 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 75, + 76 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 76, + 77 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 77, + 78 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 78, + 79 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 79, + 80 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 80, + 81 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 81, + 82 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 82, + 83 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 83, + 84 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 84, + 85 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 85, + 86 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 86, + 87 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 87, + 88 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 88, + 89 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 90, + 91 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 91, + 92 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 92, + 93 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 93, + 94 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 94, + 95 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 95, + 96 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 96, + 97 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 97, + 98 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 98, + 99 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 99, + 100 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 100, + 101 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 101, + 102 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 102, + 103 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 103, + 104 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 105, + 106 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 106, + 107 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 107, + 108 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 108, + 109 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 109, + 110 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 110, + 111 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 111, + 112 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 112, + 113 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 113, + 114 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 114, + 115 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 115, + 116 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 116, + 117 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 117, + 118 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + }, + { + "harmonic": { + "index": [ + 118, + 119 + ], + "k": 0.2758300000000001, + "req": 6.0 + } + } + ], + "relative time": 8.997733278574468e-05 + } + }, + { + "nonbonded": { + "functor potential": { + "default": [ + { + "lennardjones": { + "mixing": "LB" + } + } + ] + }, + "relative time": 0.28360199064801583, + "selfenergy": { + "dipole": false, + "monopole": false + }, + "summation_policy": "serial" + } + } + ] + } + ], + "geometry": { + "radius": 75.0, + "type": "sphere" + }, + "git revision": "3d460516 (2024-07-24)", + "groups": [ + { + "polymer": { + "compressible": false, + "index": [ + 0, + 14 + ], + "size": 15 + } + }, + { + "polymer": { + "compressible": false, + "index": [ + 15, + 29 + ], + "size": 15 + } + }, + { + "polymer": { + "compressible": false, + "index": [ + 30, + 44 + ], + "size": 15 + } + }, + { + "polymer": { + "compressible": false, + "index": [ + 45, + 59 + ], + "size": 15 + } + }, + { + "polymer": { + "compressible": false, + "index": [ + 60, + 74 + ], + "size": 15 + } + }, + { + "polymer": { + "compressible": false, + "index": [ + 75, + 89 + ], + "size": 15 + } + }, + { + "polymer": { + "compressible": false, + "index": [ + 90, + 104 + ], + "size": 15 + } + }, + { + "polymer": { + "compressible": false, + "index": [ + 105, + 119 + ], + "size": 15 + } + } + ], + "montecarlo": { + "average potential energy (kT)": 86.6550243271054, + "last move": "transrot" + }, + "moves": [ + { + "moltransrot": { + "R/Å ≈ √(⟨r²⟩/4⟨θ²⟩)": 9.77, + "acceptance": 0.78, + "dir": [ + 1.0, + 1.0, + 1.0 + ], + "dirrot": [ + 0.0, + 0.0, + 0.0 + ], + "dp": 10.0, + "dprot": 1.0, + "molecule": "polymer", + "molid": 0, + "moves": 79530, + "relative time": 0.0614, + "repeat": 8, + "stochastic": true, + "√⟨r²⟩": 4.91, + "√⟨θ²⟩": 0.251, + "√⟨θ²⟩/°": 14.4 + } + }, + { + "pivot": { + "acceptance": 0.708, + "dprot": 3.0, + "molecule": "polymer", + "moves": 1120447, + "relative time": 0.458, + "repeat": 112, + "stochastic": true, + "√⟨r_cm²⟩": 5.51 + } + }, + { + "transrot": { + "acceptance": 0.399, + "dir": [ + 1.0, + 1.0, + 1.0 + ], + "dp": 0.0, + "dprot": 0.0, + "molecule": "polymer", + "molid": 0, + "moves": 1200023, + "relative time": 0.0968, + "repeat": 120, + "stochastic": true, + "√⟨r²⟩": 2.57 + } + } + ], + "number of groups": 8, + "number of particles": 120, + "number of sweeps": 10000, + "reactionlist": null, + "relative drift": -4.3779073032648643e-14, + "simulation time": { + "in minutes": 1.2666666666666666, + "in seconds": 76 + }, + "temperature": 298.0 +} diff --git a/examples/voronota/voronota_polymers.yml b/examples/voronota/voronota_polymers.yml new file mode 100755 index 000000000..884d353f5 --- /dev/null +++ b/examples/voronota/voronota_polymers.yml @@ -0,0 +1,83 @@ +#!/usr/bin/env yason.py +random: {seed: fixed} +temperature: 298 +geometry: {type: sphere, radius: 75} +mcloop: {macro: 10, micro: 1000} + +atomlist: + - MM: {sigma: 6, dp: 10, q: 1.0, eps: 0.2} + +moleculelist: + - polymer: + structure: + - MM: [0.0, 0.0, 0.0] + - MM: [7.6, 0.0, 0.0] + - MM: [0.0, 7.6, 0.0] + - MM: [7.6, 7.6, 0.0] + - MM: [0.0, 0.0, 7.6] + - MM: [7.6, 0.0, 7.6] + - MM: [0.0, 7.6, 7.6] + - MM: [7.6, 7.6, 7.6] + - MM: [15.2, 7.6, 7.6] + - MM: [7.6, 15.2, 7.6] + - MM: [15.2, 15.2, 7.6] + - MM: [7.6, 7.6, 15.2] + - MM: [15.2, 7.6, 15.2] + - MM: [7.6, 15.2, 15.2] + - MM: [15.2, 15.2, 15.2] + bondlist: + - harmonic: {index: [0,1], k: 0.27583, req: 6} + - harmonic: {index: [1,2], k: 0.27583, req: 6} + - harmonic: {index: [2,3], k: 0.27583, req: 6} + - harmonic: {index: [3,4], k: 0.27583, req: 6} + - harmonic: {index: [4,5], k: 0.27583, req: 6} + - harmonic: {index: [5,6], k: 0.27583, req: 6} + - harmonic: {index: [6,7], k: 0.27583, req: 6} + - harmonic: {index: [7,8], k: 0.27583, req: 6} + - harmonic: {index: [8,9], k: 0.27583, req: 6} + - harmonic: {index: [9,10], k: 0.27583, req: 6} + - harmonic: {index: [10,11], k: 0.27583, req: 6} + - harmonic: {index: [11,12], k: 0.27583, req: 6} + - harmonic: {index: [12,13], k: 0.27583, req: 6} + - harmonic: {index: [13,14], k: 0.27583, req: 6} + +insertmolecules: + - polymer: {N: 8} + +energy: + - bonded: {} + - nonbonded: + default: + - lennardjones: {mixing: LB} + +moves: + - moltransrot: {molecule: polymer, dp: 10, dprot: 1, repeat: N} + - pivot: {molecule: polymer, dprot: 3, repeat: N} + - transrot: {molecule: polymer, repeat: N} + +analysis: + - sanity: {nstep: 10} + - xtcfile: {file: polymers_traj.xtc, nstep: 1} + - savestate: {file: polymers_structure.pqr} + - sasa: + file: polymers_sasa.dat + policy: molecular + molecule: polymer + radius: 2.0 + nstep: 1 + - voronoi: + file: polymers_voronoi_full.dat + mode: full + radius: 2.0 + nstep: 1 + - voronoi: + file: polymers_voronoi_interchain.dat + mode: interchain + radius: 2.0 + nstep: 1 + - voronoi: + file: polymers_voronoi_updateable.dat + mode: updateable + radius: 2.0 + nstep: 1 + diff --git a/src/analysis.h b/src/analysis.h index d658aeb32..fc0e37f94 100644 --- a/src/analysis.h +++ b/src/analysis.h @@ -1208,19 +1208,33 @@ class SavePenaltyEnergy : public Analysis */ class Voronota : public Analysis { + public: + enum class Modes + { + FULL, + INTERCHAIN, + UPDATEABLE, + INVALID + }; + private: struct Averages { - Average area; - Average area_squared; + Average area_of_contacts; + Average volume_inside_sas; + Average area_of_sas; + Average area_of_sas_per_chain; }; //!< Placeholder class for average properties - Averages average_data; //!< Stores all averages for the selected molecule - - std::unique_ptr output_stream; //!< output stream + Modes mode = Modes::FULL; //!< running mode double probe_radius; //!< radius of the probe sphere + bool use_pbc = false; //!< is the cell periodic? + Averages average_data; //!< stores all averages for the selected molecule + std::unique_ptr output_stream; //!< output stream std::string filename; //!< output file name - bool use_pbc = false; //!< Is the cell periodic? + + class Impl; + std::unique_ptr pimpl; void _to_json(json& json_output) const override; void _from_json(const json& input) override; @@ -1229,7 +1243,8 @@ class Voronota : public Analysis public: Voronota(const json& j, const Space& spc); - Voronota(double probe_radius, const Space& spc); + Voronota(Modes mode, double probe_radius, const Space& spc); + ~Voronota(); }; } // namespace Faunus::analysis diff --git a/src/voronota.cpp b/src/voronota.cpp index b06d1ef71..a1f141b83 100644 --- a/src/voronota.cpp +++ b/src/voronota.cpp @@ -6,48 +6,193 @@ namespace Faunus::analysis { -Voronota::Voronota(double probe_radius, const Faunus::Space& spc) +voronotalt::PeriodicBox get_periodic_box_from_space(const Faunus::Space& spc) +{ + const Point corner = 0.5 * spc.geometry.getLength(); + const std::vector box_corners = { + voronotalt::SimplePoint(0.0 - corner.x(), 0.0 - corner.y(), 0.0 - corner.z()), + voronotalt::SimplePoint(corner.x(), corner.y(), corner.z())}; + return voronotalt::PeriodicBox::create_periodic_box_from_corners(box_corners); +} + +Voronota::Modes get_voronota_mode_from_name(const std::string& name) +{ + if (name == "full") { + return Voronota::Modes::FULL; + } + else if (name == "interchain") { + return Voronota::Modes::INTERCHAIN; + } + else if (name == "updateable") { + return Voronota::Modes::UPDATEABLE; + } + return Voronota::Modes::INVALID; +} + +std::string get_name_from_voronota_mode(const Voronota::Modes mode) +{ + if (mode == Voronota::Modes::FULL) { + return "full"s; + } + else if (mode == Voronota::Modes::INTERCHAIN) { + return "interchain"s; + } + else if (mode == Voronota::Modes::UPDATEABLE) { + return "updateable"s; + } + return "invalid"s; +} + +class Voronota::Impl +{ + public: + Impl(const Voronota::Modes mode, const double probe_radius, const bool use_pbc) + : mode(mode) + , probe_radius(probe_radius) + , use_pbc(use_pbc) + , updateable_tessellation(false) + , num_of_chains(0) + { + } + + void reset_all(const Faunus::Space& spc) + { + periodic_box = voronotalt::PeriodicBox(); + spheres.clear(); + grouping_of_spheres.clear(); + num_of_chains = 0; + + if (use_pbc) { + periodic_box = get_periodic_box_from_space(spc); + } + + for (const auto& group : spc.groups) { + bool new_chain_active = false; + for (const auto& particle : group) { + if (!new_chain_active) { + num_of_chains++; + new_chain_active = true; + } + spheres.push_back( + voronotalt::SimpleSphere(particle.pos.x(), particle.pos.y(), particle.pos.z(), + (0.5 * particle.traits().sigma) + probe_radius)); + if (mode == Voronota::Modes::INTERCHAIN) { + grouping_of_spheres.push_back(num_of_chains); + } + } + } + + if (mode == Voronota::Modes::UPDATEABLE) { + updateable_tessellation.init(spheres, periodic_box); + } + } + + void reset_spheres(const Faunus::Space& spc) + { + if ((use_pbc != periodic_box.enabled()) || + (use_pbc && periodic_box.enabled() && + !periodic_box.equals(get_periodic_box_from_space(spc)))) { + reset_all(spc); + return; + } + + int current_num_of_chains = 0; + std::size_t sphere_index = 0; + for (const auto& group : spc.groups) { + bool new_chain_active = false; + for (const auto& particle : group) { + if (!new_chain_active) { + current_num_of_chains++; + new_chain_active = true; + } + if ((sphere_index >= spheres.size()) || + (mode == Voronota::Modes::INTERCHAIN && + (sphere_index >= grouping_of_spheres.size() || + grouping_of_spheres[sphere_index] != current_num_of_chains))) { + reset_all(spc); + return; + } + voronotalt::SimpleSphere& s = spheres[sphere_index]; + s.p.x = particle.pos.x(); + s.p.y = particle.pos.y(); + s.p.z = particle.pos.z(); + s.r = (0.5 * particle.traits().sigma) + probe_radius; + sphere_index++; + } + } + } + + Voronota::Modes mode; + double probe_radius; + bool use_pbc; + voronotalt::PeriodicBox periodic_box; + std::vector spheres; + std::vector grouping_of_spheres; + voronotalt::UpdateableRadicalTessellation updateable_tessellation; + int num_of_chains; +}; + +Voronota::Voronota(Voronota::Modes mode, double probe_radius, const Faunus::Space& spc) : Analysis(spc, "voronota") + , mode(mode) , probe_radius(probe_radius) + , use_pbc(false) { cite = "doi:10/mq8k"; + + if (mode == Voronota::Modes::INVALID) { + throw ConfigurationError("invalid running mode"); + } + auto n_pbc = spc.geometry.asSimpleGeometry()->boundary_conditions.isPeriodic().count(); - switch (n_pbc) { - case 0: - faunus_logger->debug("{}: No PBC detected", name); - use_pbc = false; - break; - case 3: - faunus_logger->debug("{}: 3D PBC detected", name); + if (n_pbc == 3) { use_pbc = true; - break; - default: + faunus_logger->debug("{}: 3D PBC detected", name); + } + else if (n_pbc == 0) { + faunus_logger->debug("{}: No PBC detected", name); + } + else { faunus_logger->warn("{}: Non-uniform PBC is currently ignored - be careful!", name); - use_pbc = false; } } Voronota::Voronota(const Faunus::json& input, const Faunus::Space& spc) - : Voronota(input.value("radius", 1.4_angstrom), spc) + : Voronota(get_voronota_mode_from_name(input.value("mode", "full"s)), + input.value("radius", 1.4_angstrom), spc) { from_json(input); } +Voronota::~Voronota() = default; + void Voronota::_from_json(const Faunus::json& input) { if (filename = input.value("file", ""s); !filename.empty()) { output_stream = IO::openCompressedOutputStream(MPI::prefix + filename, true); - *output_stream << "# step SASA\n"; + *output_stream << "# step chains spheres "; + if (mode == Voronota::Modes::FULL) { + *output_stream << "collisions relevant_collisions area_of_contacts volume_inside_sas " + "area_of_sas area_of_sas_per_chain"; + } + else if (mode == Voronota::Modes::INTERCHAIN) { + *output_stream << "collisions relevant_collisions area_of_contacts"; + } + else if (mode == Voronota::Modes::UPDATEABLE) { + *output_stream + << " area_of_contacts volume_inside_sas area_of_sas area_of_sas_per_chain"; + } + *output_stream << "\n"; } } void Voronota::_to_json(json& json_output) const { - if (!average_data.area.empty()) { - json_output = {{"⟨SASA⟩", average_data.area.avg()}, - {"⟨SASA²⟩-⟨SASA⟩²", - average_data.area_squared.avg() - std::pow(average_data.area.avg(), 2)}}; - } + json_output = {{"⟨area_of_contacts⟩", average_data.area_of_contacts.avg()}, + {"⟨volume_inside_sas⟩", average_data.volume_inside_sas.avg()}, + {"⟨area_of_sas⟩", average_data.area_of_sas.avg()}, + {"⟨area_of_sas_per_chain⟩", average_data.area_of_sas_per_chain.avg()}}; + json_output["mode"] = get_name_from_voronota_mode(mode); json_output["radius"] = probe_radius; } @@ -60,34 +205,76 @@ void Voronota::_to_disk() void Voronota::_sample() { - using voronotalt::SimplePoint; - using namespace std::views; + if (!pimpl) { + pimpl = std::make_unique(mode, probe_radius, use_pbc); + pimpl->reset_all(spc); + } + else { + pimpl->reset_spheres(spc); + } + + if (mode == Voronota::Modes::FULL) { + voronotalt::RadicalTessellation::Result result; + voronotalt::RadicalTessellation::construct_full_tessellation(pimpl->spheres, + pimpl->periodic_box, result); - // Convert single `Particle` to Voronota's `SimpleSphere` - auto to_sphere = [&](const Particle& p) -> voronotalt::SimpleSphere { - return {{p.pos.x(), p.pos.y(), p.pos.z()}, (0.5 * p.traits().sigma) + probe_radius}; - }; - const auto spheres = spc.activeParticles() | transform(to_sphere) | ranges::to_vector; + const double area_of_contacts = result.total_contacts_summary.area; + const double volume_inside_sas = result.total_cells_summary.sas_inside_volume; + const double area_of_sas = result.total_cells_summary.sas_area; + const double area_of_sas_per_chain = + (area_of_sas / static_cast(pimpl->num_of_chains)); - voronotalt::RadicalTessellation::Result result; + average_data.area_of_contacts.add(area_of_contacts); + average_data.volume_inside_sas.add(volume_inside_sas); + average_data.area_of_sas.add(area_of_sas); + average_data.area_of_sas_per_chain.add(area_of_sas_per_chain); - if (use_pbc) { - auto to_point = [](const Point& p) -> SimplePoint { return {p.x(), p.y(), p.z()}; }; - const Point corner = 0.5 * spc.geometry.getLength(); - const std::vector box_corners = {to_point(-corner), to_point(corner)}; - voronotalt::RadicalTessellation::construct_full_tessellation(spheres, box_corners, result); + if (output_stream) { + *output_stream << this->getNumberOfSteps() << " " << pimpl->num_of_chains << " " + << pimpl->spheres.size() << " " << result.total_collisions << " " + << result.total_relevant_collisions << " " << area_of_contacts << " " + << volume_inside_sas << " " << area_of_sas << " " + << area_of_sas_per_chain << "\n"; + } } - else { - voronotalt::RadicalTessellation::construct_full_tessellation(spheres, result); + else if (mode == Voronota::Modes::INTERCHAIN) { + voronotalt::RadicalTessellation::Result result; + voronotalt::RadicalTessellation::construct_full_tessellation( + pimpl->spheres, pimpl->grouping_of_spheres, pimpl->periodic_box, result); + + const double area_of_contacts = result.total_contacts_summary.area; + + average_data.area_of_contacts.add(area_of_contacts); + + if (output_stream) { + *output_stream << this->getNumberOfSteps() << " " << pimpl->num_of_chains << " " + << pimpl->spheres.size() << " " << result.total_collisions << " " + << result.total_relevant_collisions << " " << area_of_contacts << "\n"; + } } + else if (mode == Voronota::Modes::UPDATEABLE) { + pimpl->updateable_tessellation.update(pimpl->spheres); - auto areas = result.cells_summaries | transform([](auto& summary) { return summary.sas_area; }); - const auto total_area = std::accumulate(areas.begin(), areas.end(), 0.0); - average_data.area.add(total_area); - average_data.area_squared.add(total_area * total_area); + voronotalt::UpdateableRadicalTessellation::ResultSummary result = + pimpl->updateable_tessellation.result_summary(); - if (output_stream) { - *output_stream << this->getNumberOfSteps() << " " << total_area << "\n"; + const double area_of_contacts = result.total_contacts_summary.area; + const double volume_inside_sas = result.total_cells_summary.sas_inside_volume; + const double area_of_sas = result.total_cells_summary.sas_area; + const double area_of_sas_per_chain = + (area_of_sas / static_cast(pimpl->num_of_chains)); + + average_data.area_of_contacts.add(area_of_contacts); + average_data.volume_inside_sas.add(volume_inside_sas); + average_data.area_of_sas.add(area_of_sas); + average_data.area_of_sas_per_chain.add(area_of_sas_per_chain); + + if (output_stream) { + *output_stream << this->getNumberOfSteps() << " " << pimpl->num_of_chains << " " + << pimpl->spheres.size() << " " << " " << area_of_contacts << " " + << volume_inside_sas << " " << area_of_sas << " " + << area_of_sas_per_chain << "\n"; + } } } diff --git a/src/voronotalt/VERSION.txt b/src/voronotalt/VERSION.txt index c561876cf..aa130cb11 100644 --- a/src/voronotalt/VERSION.txt +++ b/src/voronotalt/VERSION.txt @@ -1 +1 @@ -Voronota-LT version 0.9.2 +Voronota-LT version 0.9.4 diff --git a/src/voronotalt/basic_types_and_functions.h b/src/voronotalt/basic_types_and_functions.h index e82d33ad6..f1615d6cc 100644 --- a/src/voronotalt/basic_types_and_functions.h +++ b/src/voronotalt/basic_types_and_functions.h @@ -37,11 +37,11 @@ struct SimplePoint Float y; Float z; - SimplePoint() : x(FLOATCONST(0.0)), y(FLOATCONST(0.0)), z(FLOATCONST(0.0)) + SimplePoint() noexcept : x(FLOATCONST(0.0)), y(FLOATCONST(0.0)), z(FLOATCONST(0.0)) { } - SimplePoint(const Float x, const Float y, const Float z) : x(x), y(y), z(z) + SimplePoint(const Float x, const Float y, const Float z) noexcept : x(x), y(y), z(z) { } }; @@ -51,11 +51,15 @@ struct SimpleSphere SimplePoint p; Float r; - SimpleSphere() : r(FLOATCONST(0.0)) + SimpleSphere() noexcept : r(FLOATCONST(0.0)) { } - SimpleSphere(const SimplePoint& p, const Float r) : p(p), r(r) + SimpleSphere(const SimplePoint& p, const Float r) noexcept : p(p), r(r) + { + } + + SimpleSphere(const Float x, const Float y, const Float z, const Float r) noexcept : p(x, y, z), r(r) { } }; @@ -67,46 +71,65 @@ struct SimpleQuaternion Float c; Float d; - SimpleQuaternion(const Float a, const Float b, const Float c, const Float d) : a(a), b(b), c(c), d(d) + SimpleQuaternion(const Float a, const Float b, const Float c, const Float d) noexcept : a(a), b(b), c(c), d(d) + { + } + + SimpleQuaternion(const Float a, const SimplePoint& p) noexcept : a(a), b(p.x), c(p.y), d(p.z) + { + } +}; + +struct ValuedID +{ + Float value; + UnsignedInt index; + + ValuedID() noexcept : value(FLOATCONST(0.0)), index(0) { } - SimpleQuaternion(const Float a, const SimplePoint& p) : a(a), b(p.x), c(p.y), d(p.z) + ValuedID(const Float value, const UnsignedInt index) noexcept : value(value), index(index) { } + + bool operator<(const ValuedID& cid) const noexcept + { + return (valueb); } -inline bool less_or_equal(const Float a, const Float b) +inline bool less_or_equal(const Float a, const Float b) noexcept { return (less(a, b) || equal(a, b)); } -inline bool greater_or_equal(const Float a, const Float b) +inline bool greater_or_equal(const Float a, const Float b) noexcept { return (greater(a, b) || equal(a, b)); } -inline void set_close_to_equal(Float& a, const Float b) +inline void set_close_to_equal(Float& a, const Float b) noexcept { if(equal(a, b)) { @@ -114,7 +137,7 @@ inline void set_close_to_equal(Float& a, const Float b) } } -inline Float squared_distance_from_point_to_point(const SimplePoint& a, const SimplePoint& b) +inline Float squared_distance_from_point_to_point(const SimplePoint& a, const SimplePoint& b) noexcept { const Float dx=(a.x-b.x); const Float dy=(a.y-b.y); @@ -122,57 +145,62 @@ inline Float squared_distance_from_point_to_point(const SimplePoint& a, const Si return (dx*dx+dy*dy+dz*dz); } -inline Float distance_from_point_to_point(const SimplePoint& a, const SimplePoint& b) +inline Float distance_from_point_to_point(const SimplePoint& a, const SimplePoint& b) noexcept { return sqrt(squared_distance_from_point_to_point(a, b)); } -inline Float squared_point_module(const SimplePoint& a) +inline bool point_equals_point(const SimplePoint& a, const SimplePoint& b) noexcept +{ + return (equal(a.x, b.x) && equal(a.y, b.y) && equal(a.z, b.z)); +} + +inline Float squared_point_module(const SimplePoint& a) noexcept { return (a.x*a.x+a.y*a.y+a.z*a.z); } -inline Float point_module(const SimplePoint& a) +inline Float point_module(const SimplePoint& a) noexcept { return sqrt(squared_point_module(a)); } -inline Float dot_product(const SimplePoint& a, const SimplePoint& b) +inline Float dot_product(const SimplePoint& a, const SimplePoint& b) noexcept { return (a.x*b.x+a.y*b.y+a.z*b.z); } -inline SimplePoint cross_product(const SimplePoint& a, const SimplePoint& b) +inline SimplePoint cross_product(const SimplePoint& a, const SimplePoint& b) noexcept { return SimplePoint(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x); } -inline SimplePoint point_and_number_product(const SimplePoint& a, const Float k) +inline SimplePoint point_and_number_product(const SimplePoint& a, const Float k) noexcept { return SimplePoint(a.x*k, a.y*k, a.z*k); } -inline SimplePoint unit_point(const SimplePoint& a) +inline SimplePoint unit_point(const SimplePoint& a) noexcept { return ((equal(squared_point_module(a), FLOATCONST(1.0))) ? a : point_and_number_product(a, FLOATCONST(1.0)/point_module(a))); } -inline SimplePoint sum_of_points(const SimplePoint& a, const SimplePoint& b) +inline SimplePoint sum_of_points(const SimplePoint& a, const SimplePoint& b) noexcept { return SimplePoint(a.x+b.x, a.y+b.y, a.z+b.z); } -inline SimplePoint sub_of_points(const SimplePoint& a, const SimplePoint& b) +inline SimplePoint sub_of_points(const SimplePoint& a, const SimplePoint& b) noexcept { return SimplePoint(a.x-b.x, a.y-b.y, a.z-b.z); } -inline Float signed_distance_from_point_to_plane(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& x) +inline Float signed_distance_from_point_to_plane(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& x) noexcept { return dot_product(unit_point(plane_normal), sub_of_points(x, plane_point)); } -inline int halfspace_of_point(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& x) +inline int halfspace_of_point(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& x) noexcept { const Float sd=signed_distance_from_point_to_plane(plane_point, plane_normal, x); if(greater(sd, FLOATCONST(0.0))) @@ -186,7 +214,7 @@ inline int halfspace_of_point(const SimplePoint& plane_point, const SimplePoint& return 0; } -inline SimplePoint intersection_of_plane_and_segment(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& a, const SimplePoint& b) +inline SimplePoint intersection_of_plane_and_segment(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& a, const SimplePoint& b) noexcept { const Float da=signed_distance_from_point_to_plane(plane_point, plane_normal, a); const Float db=signed_distance_from_point_to_plane(plane_point, plane_normal, b); @@ -201,12 +229,12 @@ inline SimplePoint intersection_of_plane_and_segment(const SimplePoint& plane_po } } -inline Float triangle_area(const SimplePoint& a, const SimplePoint& b, const SimplePoint& c) +inline Float triangle_area(const SimplePoint& a, const SimplePoint& b, const SimplePoint& c) noexcept { return (point_module(cross_product(sub_of_points(b, a), sub_of_points(c, a)))/FLOATCONST(2.0)); } -inline Float min_angle(const SimplePoint& o, const SimplePoint& a, const SimplePoint& b) +inline Float min_angle(const SimplePoint& o, const SimplePoint& a, const SimplePoint& b) noexcept { Float cos_val=dot_product(unit_point(sub_of_points(a, o)), unit_point(sub_of_points(b, o))); if(cos_valFLOATCONST(0.0)) @@ -326,7 +354,7 @@ inline bool intersect_segment_with_circle(const SimpleSphere& circle, const Simp return false; } -inline Float min_dihedral_angle(const SimplePoint& o, const SimplePoint& a, const SimplePoint& b1, const SimplePoint& b2) +inline Float min_dihedral_angle(const SimplePoint& o, const SimplePoint& a, const SimplePoint& b1, const SimplePoint& b2) noexcept { const SimplePoint oa=unit_point(sub_of_points(a, o)); const SimplePoint d1=sub_of_points(b1, sum_of_points(o, point_and_number_product(oa, dot_product(oa, sub_of_points(b1, o))))); @@ -335,7 +363,7 @@ inline Float min_dihedral_angle(const SimplePoint& o, const SimplePoint& a, cons return std::acos(std::max(FLOATCONST(-1.0), std::min(cos_val, FLOATCONST(1.0)))); } -inline SimpleQuaternion quaternion_product(const SimpleQuaternion& q1, const SimpleQuaternion& q2) +inline SimpleQuaternion quaternion_product(const SimpleQuaternion& q1, const SimpleQuaternion& q2) noexcept { return SimpleQuaternion( q1.a*q2.a - q1.b*q2.b - q1.c*q2.c - q1.d*q2.d, @@ -344,12 +372,12 @@ inline SimpleQuaternion quaternion_product(const SimpleQuaternion& q1, const Sim q1.a*q2.d + q1.b*q2.c - q1.c*q2.b + q1.d*q2.a); } -inline SimpleQuaternion inverted_quaternion(const SimpleQuaternion& q) +inline SimpleQuaternion inverted_quaternion(const SimpleQuaternion& q) noexcept { return SimpleQuaternion(q.a, FLOATCONST(0.0)-q.b, FLOATCONST(0.0)-q.c, FLOATCONST(0.0)-q.d); } -inline SimplePoint rotate_point_around_axis(const SimplePoint& axis, const Float angle, const SimplePoint& p) +inline SimplePoint rotate_point_around_axis(const SimplePoint& axis, const Float angle, const SimplePoint& p) noexcept { if(squared_point_module(axis)>0) { diff --git a/src/voronotalt/conversion_to_input.h b/src/voronotalt/conversion_to_input.h index 473530f6c..bbf6eeaca 100644 --- a/src/voronotalt/conversion_to_input.h +++ b/src/voronotalt/conversion_to_input.h @@ -9,13 +9,13 @@ namespace voronotalt { template -inline SimpleSphere get_sphere_from_ball(const Ball& ball, const Float probe) +inline SimpleSphere get_sphere_from_ball(const Ball& ball, const Float probe) noexcept { return SimpleSphere(SimplePoint(static_cast(ball.x), static_cast(ball.y), static_cast(ball.z)), static_cast(ball.r)+probe); } template -inline void fill_sphere_from_ball(const Ball& ball, const Float probe, SimpleSphere& sphere) +inline void fill_sphere_from_ball(const Ball& ball, const Float probe, SimpleSphere& sphere) noexcept { sphere.p.x=static_cast(ball.x); sphere.p.y=static_cast(ball.y); @@ -24,7 +24,7 @@ inline void fill_sphere_from_ball(const Ball& ball, const Float probe, SimpleSph } template -inline std::vector get_spheres_from_balls(const BallsContainer& balls, const Float probe) +inline std::vector get_spheres_from_balls(const BallsContainer& balls, const Float probe) noexcept { std::vector result; result.reserve(balls.size()); @@ -36,7 +36,7 @@ inline std::vector get_spheres_from_balls(const BallsContainer& ba } template -inline void fill_spheres_from_balls(const BallsContainer& balls, const Float probe, std::vector& spheres) +inline void fill_spheres_from_balls(const BallsContainer& balls, const Float probe, std::vector& spheres) noexcept { spheres.resize(balls.size()); std::size_t i=0; @@ -48,13 +48,13 @@ inline void fill_spheres_from_balls(const BallsContainer& balls, const Float pro } template -inline SimplePoint get_simple_point_from_point(const Point& point) +inline SimplePoint get_simple_point_from_point(const Point& point) noexcept { return SimplePoint(static_cast(point.x), static_cast(point.y), static_cast(point.z)); } template -inline std::vector get_simple_points_from_points(const PointsContainer& points) +inline std::vector get_simple_points_from_points(const PointsContainer& points) noexcept { std::vector result; result.reserve(points.size()); diff --git a/src/voronotalt/parallelization_configuration.h b/src/voronotalt/parallelization_configuration.h new file mode 100644 index 000000000..481b372b1 --- /dev/null +++ b/src/voronotalt/parallelization_configuration.h @@ -0,0 +1,8 @@ +#ifndef VORONOTALT_PARALLELIZATION_CONFIGURATION_H_ +#define VORONOTALT_PARALLELIZATION_CONFIGURATION_H_ + +#ifdef _OPENMP +#define VORONOTALT_OPENMP +#endif + +#endif /* VORONOTALT_PARALLELIZATION_CONFIGURATION_H_ */ diff --git a/src/voronotalt/periodic_box.h b/src/voronotalt/periodic_box.h new file mode 100644 index 000000000..583c9fcec --- /dev/null +++ b/src/voronotalt/periodic_box.h @@ -0,0 +1,94 @@ +#ifndef VORONOTALT_PERIODIC_BOX_H_ +#define VORONOTALT_PERIODIC_BOX_H_ + +#include + +#include "basic_types_and_functions.h" + +namespace voronotalt +{ + +class PeriodicBox +{ +public: + PeriodicBox() noexcept : enabled_(false) + { + } + + bool enabled() const noexcept + { + return enabled_; + } + + bool equals(const PeriodicBox& pb) const noexcept + { + return (enabled_==pb.enabled_ && point_equals_point(shift_direction_a_, pb.shift_direction_a_) && point_equals_point(shift_direction_b_, pb.shift_direction_b_) && point_equals_point(shift_direction_c_, pb.shift_direction_c_)); + } + + const SimpleSphere shift_by_weighted_directions(const SimpleSphere& o, const Float wa, const Float wb, const Float wc) const + { + return SimpleSphere( + o.p.x+shift_direction_a_.x*wa+shift_direction_b_.x*wb+shift_direction_c_.x*wc, + o.p.y+shift_direction_a_.y*wa+shift_direction_b_.y*wb+shift_direction_c_.y*wc, + o.p.z+shift_direction_a_.z*wa+shift_direction_b_.z*wb+shift_direction_c_.z*wc, + o.r); + } + + static PeriodicBox create_periodic_box_from_corners(const std::vector& periodic_box_corners) noexcept + { + PeriodicBox pb; + if(periodic_box_corners.size()>=2) + { + pb.enabled_=true; + SimplePoint corner_min=periodic_box_corners[0]; + SimplePoint corner_max=periodic_box_corners[0]; + for(UnsignedInt i=1;i& shift_directions) noexcept + { + PeriodicBox pb; + if(shift_directions.size()==3) + { + pb.enabled_=true; + pb.shift_direction_a_=shift_directions[0]; + pb.shift_direction_b_=shift_directions[1]; + pb.shift_direction_c_=shift_directions[2]; + } + return pb; + } + + static PeriodicBox create_periodic_box_from_shift_directions_or_from_corners(const std::vector& shift_directions, const std::vector& periodic_box_corners) noexcept + { + PeriodicBox pb=create_periodic_box_from_shift_directions(shift_directions); + if(!pb.enabled()) + { + pb=create_periodic_box_from_corners(periodic_box_corners); + } + return pb; + } + +private: + bool enabled_; + SimplePoint shift_direction_a_; + SimplePoint shift_direction_b_; + SimplePoint shift_direction_c_; +}; + +} + +#endif /* VORONOTALT_PERIODIC_BOX_H_ */ diff --git a/src/voronotalt/preparation_for_tessellation.h b/src/voronotalt/preparation_for_tessellation.h deleted file mode 100644 index 31a6a964f..000000000 --- a/src/voronotalt/preparation_for_tessellation.h +++ /dev/null @@ -1,193 +0,0 @@ -#ifndef VORONOTALT_PREPARATION_FOR_TESSELLATION_H_ -#define VORONOTALT_PREPARATION_FOR_TESSELLATION_H_ - -#include - -#include "spheres_searcher.h" -#include "time_recorder.h" - -namespace voronotalt -{ - -class PreparationForTessellation -{ -public: - struct Result - { - std::vector spheres_with_periodic_instances; - std::vector periodic_links_of_spheres; - std::vector all_exclusion_statuses; - std::vector< std::vector > all_colliding_ids; - std::vector< std::pair > relevant_collision_ids; - UnsignedInt total_input_spheres; - UnsignedInt total_spheres; - UnsignedInt total_collisions; - - Result() : total_input_spheres(0), total_spheres(0), total_collisions(0) - { - } - }; - - static void prepare_for_tessellation( - const std::vector& spheres, - const std::vector& grouping_of_spheres, - Result& result, - TimeRecorder& time_recorder) - { - prepare_for_tessellation(spheres, grouping_of_spheres, std::vector(), result, time_recorder); - } - - static void prepare_for_tessellation( - const std::vector& input_spheres, - const std::vector& grouping_of_spheres, - const std::vector& periodic_box_corners, - Result& result, - TimeRecorder& time_recorder) - { - time_recorder.reset(); - - result=Result(); - - result.total_input_spheres=input_spheres.size(); - - if(periodic_box_corners.size()>=2) - { - SimplePoint corner_a=periodic_box_corners[0]; - SimplePoint corner_b=periodic_box_corners[0]; - for(UnsignedInt i=1;i(sx)); - const bool useful_x=((sx<0 && (corner_a.x-m.p.x)<=buffer_distance) || (sx>0 && (m.p.x-corner_b.x)<=buffer_distance)); - for(int sy=-1;sy<=1;sy++) - { - m.p.y=o.p.y+(shift.y*static_cast(sy)); - const bool useful_y=((sy<0 && (corner_a.y-m.p.y)<=buffer_distance) || (sy>0 && (m.p.y-corner_b.y)<=buffer_distance)); - for(int sz=-1;sz<=1;sz++) - { - if(sx!=0 || sy!=0 || sz!=0) - { - m.p.z=o.p.z+(shift.z*static_cast(sz)); - const bool useful_z=((sz<0 && (corner_a.z-m.p.z)<=buffer_distance) || (sz>0 && (m.p.z-corner_b.z)<=buffer_distance)); - if(useful_x || useful_y || useful_z) - { - result.spheres_with_periodic_instances.push_back(m); - result.periodic_links_of_spheres.push_back(i); - } - } - } - } - } - } - } - - const std::vector& spheres=(result.spheres_with_periodic_instances.empty() ? input_spheres : result.spheres_with_periodic_instances); - - const UnsignedInt N=spheres.size(); - result.total_spheres=N; - - SpheresSearcher spheres_searcher(spheres); - - time_recorder.record_elapsed_miliseconds_and_reset("init spheres searcher"); - - result.all_exclusion_statuses.resize(N, 0); - - result.all_colliding_ids.resize(N); - for(UnsignedInt i=0;i colliding_ids; - colliding_ids.reserve(100); - - std::vector< std::pair > distances_of_colliding_ids; - distances_of_colliding_ids.reserve(100); - - #pragma omp for - for(UnsignedInt i=0;i=grouping_of_spheres.size() || id_b>=grouping_of_spheres.size() || grouping_of_spheres[id_a]!=grouping_of_spheres[id_b]) - { - if(result.periodic_links_of_spheres.empty() || id_a>=result.periodic_links_of_spheres.size() || id_b>=result.periodic_links_of_spheres.size() || result.periodic_links_of_spheres[id_a]==id_a || result.periodic_links_of_spheres[id_b]==id_b) - { - result.relevant_collision_ids.push_back(std::pair(id_a, id_b)); - } - } - } - } - } - } - - time_recorder.record_elapsed_miliseconds_and_reset("collect relevant collision indices"); - } -}; - -} - -#endif /* VORONOTALT_PREPARATION_FOR_TESSELLATION_H_ */ diff --git a/src/voronotalt/radical_tessellation.h b/src/voronotalt/radical_tessellation.h index 909147d7a..0941caeaf 100644 --- a/src/voronotalt/radical_tessellation.h +++ b/src/voronotalt/radical_tessellation.h @@ -3,7 +3,7 @@ #include -#include "preparation_for_tessellation.h" +#include "spheres_container.h" #include "radical_tessellation_contact_construction.h" #include "time_recorder.h" @@ -26,7 +26,7 @@ class RadicalTessellation UnsignedInt id_a; UnsignedInt id_b; - ContactDescriptorSummary() : + ContactDescriptorSummary() noexcept : area(FLOATCONST(0.0)), arc_length(FLOATCONST(0.0)), solid_angle_a(FLOATCONST(0.0)), @@ -40,7 +40,7 @@ class RadicalTessellation { } - void set(const RadicalTessellationContactConstruction::ContactDescriptor& cd) + void set(const RadicalTessellationContactConstruction::ContactDescriptor& cd) noexcept { if(cd.area>FLOATCONST(0.0)) { @@ -57,7 +57,7 @@ class RadicalTessellation } } - void ensure_ids_ordered() + void ensure_ids_ordered() noexcept { if(id_a>id_b) { @@ -68,6 +68,19 @@ class RadicalTessellation } }; + struct ContactDescriptorSummaryAdjunct + { + ContactDescriptorSummaryAdjunct() noexcept + { + } + + explicit ContactDescriptorSummaryAdjunct(const UnsignedInt levels) noexcept : level_areas(levels, FLOATCONST(0.0)) + { + } + + std::vector level_areas; + }; + struct CellContactDescriptorsSummary { Float area; @@ -82,7 +95,7 @@ class RadicalTessellation UnsignedInt count; int stage; - CellContactDescriptorsSummary() : + CellContactDescriptorsSummary() noexcept : area(FLOATCONST(0.0)), arc_length(FLOATCONST(0.0)), explained_solid_angle_positive(FLOATCONST(0.0)), @@ -97,7 +110,7 @@ class RadicalTessellation { } - void add(const ContactDescriptorSummary& cds) + void add(const ContactDescriptorSummary& cds) noexcept { if(cds.area>FLOATCONST(0.0) && (cds.id_a==id || cds.id_b==id)) { @@ -112,7 +125,7 @@ class RadicalTessellation } } - void add(const UnsignedInt new_id, const ContactDescriptorSummary& cds) + void add(const UnsignedInt new_id, const ContactDescriptorSummary& cds) noexcept { if(cds.area>FLOATCONST(0.0)) { @@ -124,7 +137,7 @@ class RadicalTessellation } } - void compute_sas(const Float r) + void compute_sas(const Float r) noexcept { if(stage==1) { @@ -155,7 +168,7 @@ class RadicalTessellation } } - void compute_sas_detached(const UnsignedInt new_id, const Float r) + void compute_sas_detached(const UnsignedInt new_id, const Float r) noexcept { if(stage==0) { @@ -174,7 +187,7 @@ class RadicalTessellation Float distance; UnsignedInt count; - TotalContactDescriptorsSummary() : + TotalContactDescriptorsSummary() noexcept : area(FLOATCONST(0.0)), arc_length(FLOATCONST(0.0)), distance(FLOATCONST(-1.0)), @@ -182,7 +195,7 @@ class RadicalTessellation { } - void add(const ContactDescriptorSummary& cds) + void add(const ContactDescriptorSummary& cds) noexcept { if(cds.area>FLOATCONST(0.0)) { @@ -200,14 +213,14 @@ class RadicalTessellation Float sas_inside_volume; UnsignedInt count; - TotalCellContactDescriptorsSummary() : + TotalCellContactDescriptorsSummary() noexcept : sas_area(FLOATCONST(0.0)), sas_inside_volume(FLOATCONST(0.0)), count(0) { } - void add(const CellContactDescriptorsSummary& ccds) + void add(const CellContactDescriptorsSummary& ccds) noexcept { if(ccds.stage==2) { @@ -224,14 +237,29 @@ class RadicalTessellation UnsignedInt total_collisions; UnsignedInt total_relevant_collisions; std::vector contacts_summaries; + std::vector adjuncts_for_contacts_summaries; TotalContactDescriptorsSummary total_contacts_summary; std::vector cells_summaries; TotalCellContactDescriptorsSummary total_cells_summary; - std::vector contacts_with_redundancy_summaries_in_periodic_box; - std::vector contacts_with_redundancy_canonical_ids_in_periodic_box; + std::vector contacts_summaries_with_redundancy_in_periodic_box; + std::vector contacts_canonical_ids_with_redundancy_in_periodic_box; + + Result() noexcept : total_spheres(0), total_collisions(0), total_relevant_collisions(0) + { + } - Result() : total_spheres(0), total_collisions(0), total_relevant_collisions(0) + void clear() noexcept { + total_spheres=0; + total_collisions=0; + total_relevant_collisions=0; + contacts_summaries.clear(); + adjuncts_for_contacts_summaries.clear(); + total_contacts_summary=TotalContactDescriptorsSummary(); + cells_summaries.clear(); + total_cells_summary=TotalCellContactDescriptorsSummary(); + contacts_summaries_with_redundancy_in_periodic_box.clear(); + contacts_canonical_ids_with_redundancy_in_periodic_box.clear(); } }; @@ -246,49 +274,99 @@ class RadicalTessellation std::vector grouped_contacts_summaries; std::vector grouped_cells_representative_ids; std::vector grouped_cells_summaries; + + void clear() noexcept + { + grouped_contacts_representative_ids.clear(); + grouped_contacts_summaries.clear(); + grouped_cells_representative_ids.clear(); + grouped_cells_summaries.clear(); + } }; static void construct_full_tessellation( const std::vector& input_spheres, - Result& result) + Result& result) noexcept + { + TimeRecorder time_recorder; + SpheresContainer spheres_container; + spheres_container.init(input_spheres, time_recorder); + ResultGraphics result_graphics; + construct_full_tessellation(spheres_container, std::vector(), std::vector(), false, true, FLOATCONST(0.0), std::vector(), result, result_graphics, time_recorder); + } + + static void construct_full_tessellation( + const std::vector& input_spheres, + const std::vector& grouping_of_spheres, + Result& result) noexcept { TimeRecorder time_recorder; + SpheresContainer spheres_container; + spheres_container.init(input_spheres, time_recorder); ResultGraphics result_graphics; - construct_full_tessellation(input_spheres, std::vector(), std::vector(), false, true, result, result_graphics, time_recorder); + construct_full_tessellation(spheres_container, std::vector(), grouping_of_spheres, false, grouping_of_spheres.empty(), FLOATCONST(0.0), std::vector(), result, result_graphics, time_recorder); } static void construct_full_tessellation( const std::vector& input_spheres, - const std::vector& periodic_box_corners, - Result& result) + const PeriodicBox& periodic_box, + Result& result) noexcept { TimeRecorder time_recorder; + SpheresContainer spheres_container; + spheres_container.init(input_spheres, periodic_box, time_recorder); ResultGraphics result_graphics; - construct_full_tessellation(input_spheres, std::vector(), periodic_box_corners, false, true, result, result_graphics, time_recorder); + construct_full_tessellation(spheres_container, std::vector(), std::vector(), false, true, FLOATCONST(0.0), std::vector(), result, result_graphics, time_recorder); } static void construct_full_tessellation( const std::vector& input_spheres, const std::vector& grouping_of_spheres, - const std::vector& periodic_box_corners, + const PeriodicBox& periodic_box, + Result& result) noexcept + { + TimeRecorder time_recorder; + SpheresContainer spheres_container; + spheres_container.init(input_spheres, periodic_box, time_recorder); + ResultGraphics result_graphics; + construct_full_tessellation(spheres_container, std::vector(), grouping_of_spheres, false, grouping_of_spheres.empty(), FLOATCONST(0.0), std::vector(), result, result_graphics, time_recorder); + } + + static void construct_full_tessellation( + const SpheresContainer& spheres_container, + const std::vector& grouping_of_spheres, const bool with_graphics, const bool summarize_cells, Result& result, ResultGraphics& result_graphics, - TimeRecorder& time_recorder) + TimeRecorder& time_recorder) noexcept { - PreparationForTessellation::Result preparation_result; - PreparationForTessellation::prepare_for_tessellation(input_spheres, grouping_of_spheres, periodic_box_corners, preparation_result, time_recorder); - - const std::vector& spheres=(preparation_result.spheres_with_periodic_instances.empty() ? input_spheres : preparation_result.spheres_with_periodic_instances); + construct_full_tessellation(spheres_container, std::vector(), grouping_of_spheres, with_graphics, summarize_cells, FLOATCONST(0.0), std::vector(), result, result_graphics, time_recorder); + } + static void construct_full_tessellation( + const SpheresContainer& spheres_container, + const std::vector& involvement_of_spheres, + const std::vector& grouping_of_spheres, + const bool with_graphics, + const bool summarize_cells, + const Float max_circle_radius_restriction, + const std::vector& adjunct_max_circle_radius_restrictions, + Result& result, + ResultGraphics& result_graphics, + TimeRecorder& time_recorder) noexcept + { time_recorder.reset(); - result=Result(); + result.clear(); + + SpheresContainer::ResultOfPreparationForTessellation preparation_result; + spheres_container.prepare_for_tessellation(involvement_of_spheres, grouping_of_spheres, preparation_result, time_recorder); + result_graphics=ResultGraphics(); - result.total_spheres=preparation_result.total_input_spheres; - result.total_collisions=preparation_result.total_collisions; + result.total_spheres=spheres_container.input_spheres().size(); + result.total_collisions=spheres_container.total_collisions(); result.total_relevant_collisions=preparation_result.relevant_collision_ids.size(); std::vector possible_contacts_summaries(preparation_result.relevant_collision_ids.size()); @@ -301,17 +379,28 @@ class RadicalTessellation time_recorder.record_elapsed_miliseconds_and_reset("allocate possible contact summaries"); - #pragma omp parallel +#ifdef VORONOTALT_OPENMP +#pragma omp parallel +#endif { RadicalTessellationContactConstruction::ContactDescriptor cd; cd.contour.reserve(12); - #pragma omp for +#ifdef VORONOTALT_OPENMP +#pragma omp for +#endif for(UnsignedInt i=0;iFLOATCONST(0.0)) + { + ContactDescriptorSummaryAdjunct& cdsa=result.adjuncts_for_contacts_summaries[i]; + Float prev_circle_radius_restriction=0.0; + for(UnsignedInt j=0;jFLOATCONST(0.0) ? std::min(adjunct_max_circle_radius_restrictions[j], max_circle_radius_restriction) : adjunct_max_circle_radius_restrictions[j]); + if(j==0 || (circle_radius_restriction>=prev_circle_radius_restriction) + || (circle_radius_restrictionFLOATCONST(0.0))) + { + if(RadicalTessellationContactConstruction::construct_contact_descriptor( + spheres_container.populated_spheres(), + spheres_container.all_exclusion_statuses(), + cds.id_a, + cds.id_b, + spheres_container.all_colliding_ids()[cds.id_a], + circle_radius_restriction, + cd)) + { + cdsa.level_areas[j]=cd.area; + } + } + prev_circle_radius_restriction=circle_radius_restriction; + } + } } } + } - time_recorder.record_elapsed_miliseconds_and_reset("accumulate cell summaries"); - - for(UnsignedInt i=0;i > map_of_spheres_to_boundary_contacts(result.total_spheres); for(UnsignedInt i=0;i=result.total_spheres || cds.id_b>=result.total_spheres) && cds.id_a=result.total_spheres || cds.id_b>=result.total_spheres) { - map_of_spheres_to_boundary_contacts[preparation_result.periodic_links_of_spheres[cds.id_a]].push_back(i); - map_of_spheres_to_boundary_contacts[preparation_result.periodic_links_of_spheres[cds.id_b]].push_back(i); + map_of_spheres_to_boundary_contacts[cds.id_a%result.total_spheres].push_back(i); + map_of_spheres_to_boundary_contacts[cds.id_b%result.total_spheres].push_back(i); } } - result.contacts_with_redundancy_canonical_ids_in_periodic_box.resize(result.contacts_summaries.size()); + result.contacts_canonical_ids_with_redundancy_in_periodic_box.resize(result.contacts_summaries.size()); UnsignedInt count_of_redundant_contacts_in_periodic_box=0; for(UnsignedInt i=0;i=result.total_spheres || cds.id_b>=result.total_spheres) - && cds.id_a=result.total_spheres || cds.id_b>=result.total_spheres) { - const UnsignedInt sphere_id_a=preparation_result.periodic_links_of_spheres[cds.id_a]; - const UnsignedInt sphere_id_b=preparation_result.periodic_links_of_spheres[cds.id_b]; + const UnsignedInt sphere_id_a=(cds.id_a%result.total_spheres); + const UnsignedInt sphere_id_b=(cds.id_b%result.total_spheres); const std::vector& candidate_ids_a=map_of_spheres_to_boundary_contacts[sphere_id_a]; const std::vector& candidate_ids_b=map_of_spheres_to_boundary_contacts[sphere_id_b]; const std::vector& candidate_ids=(candidate_ids_a.size()<=candidate_ids_b.size() ? candidate_ids_a : candidate_ids_b); @@ -447,21 +545,17 @@ class RadicalTessellation { const UnsignedInt candidate_id=candidate_ids[j]; const ContactDescriptorSummary& candidate_cds=result.contacts_summaries[candidate_id]; - if(candidate_cds.id_a0) { - result.contacts_with_redundancy_summaries_in_periodic_box.swap(result.contacts_summaries); - result.contacts_summaries.reserve(result.contacts_with_redundancy_summaries_in_periodic_box.size()+1-count_of_redundant_contacts_in_periodic_box); - for(UnsignedInt i=0;i=result.contacts_canonical_ids_with_redundancy_in_periodic_box.size() || result.contacts_canonical_ids_with_redundancy_in_periodic_box[i]==i) { - cds.id_b=preparation_result.periodic_links_of_spheres[cds.id_b]; - } - cds.ensure_ids_ordered(); - if(i>=result.contacts_with_redundancy_canonical_ids_in_periodic_box.size() || result.contacts_with_redundancy_canonical_ids_in_periodic_box[i]==i) - { - result.contacts_summaries.push_back(cds); + result.contacts_summaries.push_back(result.contacts_summaries_with_redundancy_in_periodic_box[i]); + ContactDescriptorSummary& cds=result.contacts_summaries.back(); + cds.id_a=(cds.id_a%result.total_spheres); + cds.id_b=(cds.id_b%result.total_spheres); + cds.ensure_ids_ordered(); } } } @@ -502,12 +590,52 @@ class RadicalTessellation } time_recorder.record_elapsed_miliseconds_and_reset("accumulate total contacts summary"); + + if(summarize_cells && grouping_of_spheres.empty()) + { + const std::vector& all_contacts_summaries=(result.contacts_summaries_with_redundancy_in_periodic_box.empty() ? result.contacts_summaries : result.contacts_summaries_with_redundancy_in_periodic_box); + + result.cells_summaries.resize(result.total_spheres); + + for(UnsignedInt i=0;i& grouping_of_spheres, - GroupedResult& grouped_result) + GroupedResult& grouped_result) noexcept { TimeRecorder time_recorder; return group_results(full_result, grouping_of_spheres, grouped_result, time_recorder); @@ -517,11 +645,11 @@ class RadicalTessellation const Result& full_result, const std::vector& grouping_of_spheres, GroupedResult& grouped_result, - TimeRecorder& time_recorder) + TimeRecorder& time_recorder) noexcept { time_recorder.reset(); - grouped_result=GroupedResult(); + grouped_result.clear(); if(!grouping_of_spheres.empty() && grouping_of_spheres.size()==full_result.total_spheres) { diff --git a/src/voronotalt/radical_tessellation_contact_construction.h b/src/voronotalt/radical_tessellation_contact_construction.h index 911829854..1e1993ba6 100644 --- a/src/voronotalt/radical_tessellation_contact_construction.h +++ b/src/voronotalt/radical_tessellation_contact_construction.h @@ -20,7 +20,7 @@ class RadicalTessellationContactConstruction int indicator; - ContourPoint(const SimplePoint& p, const UnsignedInt left_id, const UnsignedInt right_id) : + ContourPoint(const SimplePoint& p, const UnsignedInt left_id, const UnsignedInt right_id) noexcept : p(p), angle(FLOATCONST(0.0)), left_id(left_id), @@ -49,7 +49,7 @@ class RadicalTessellationContactConstruction UnsignedInt id_a; UnsignedInt id_b; - ContactDescriptor() : + ContactDescriptor() noexcept : sum_of_arc_angles(FLOATCONST(0.0)), area(FLOATCONST(0.0)), solid_angle_a(FLOATCONST(0.0)), @@ -63,7 +63,7 @@ class RadicalTessellationContactConstruction { } - void clear() + void clear() noexcept { id_a=0; id_b=0; @@ -85,11 +85,11 @@ class RadicalTessellationContactConstruction SimplePoint barycenter; SimplePoint plane_normal; - ContactDescriptorGraphics() + ContactDescriptorGraphics() noexcept { } - void clear() + void clear() noexcept { outer_points.clear(); } @@ -100,8 +100,9 @@ class RadicalTessellationContactConstruction const std::vector& spheres_exclusion_statuses, const UnsignedInt a_id, const UnsignedInt b_id, - const std::vector& a_neighbor_collisions, - ContactDescriptor& result_contact_descriptor) + const std::vector& a_neighbor_collisions, + const Float max_circle_radius_restriction, + ContactDescriptor& result_contact_descriptor) noexcept { result_contact_descriptor.clear(); if(a_idFLOATCONST(0.0)) + { + result_contact_descriptor.intersection_circle_sphere.r=std::min(result_contact_descriptor.intersection_circle_sphere.r, max_circle_radius_restriction); + } if(result_contact_descriptor.intersection_circle_sphere.r>FLOATCONST(0.0)) { bool discarded=false; bool contour_initialized=false; - bool nostop=true; { + bool nostop=true; for(UnsignedInt i=0;i=spheres_exclusion_statuses.size() || spheres_exclusion_statuses[neighbor_id]==0)) { const SimpleSphere& c=spheres[neighbor_id]; @@ -219,7 +224,7 @@ class RadicalTessellationContactConstruction return (result_contact_descriptor.area>FLOATCONST(0.0)); } - static bool construct_contact_descriptor_graphics(const ContactDescriptor& contact_descriptor, const Float length_step, ContactDescriptorGraphics& result_contact_descriptor_graphics) + static bool construct_contact_descriptor_graphics(const ContactDescriptor& contact_descriptor, const Float length_step, ContactDescriptorGraphics& result_contact_descriptor_graphics) noexcept { result_contact_descriptor_graphics.clear(); if(contact_descriptor.area>FLOATCONST(0.0)) @@ -291,7 +296,7 @@ class RadicalTessellationContactConstruction const UnsignedInt a_id, const SimpleSphere& base, const SimplePoint& axis, - Contour& result) + Contour& result) noexcept { const SimplePoint first_point=point_and_number_product(any_normal_of_vector(axis), base.r*FLOATCONST(1.19)); const Float angle_step=PIVALUE/FLOATCONST(3.0); @@ -306,7 +311,7 @@ class RadicalTessellationContactConstruction const SimplePoint& ac_plane_center, const SimplePoint& ac_plane_normal, const UnsignedInt c_id, - Contour& contour) + Contour& contour) noexcept { const UnsignedInt outsiders_count=mark_contour(ac_plane_center, ac_plane_normal, c_id, contour); if(outsiders_count>0) @@ -328,7 +333,7 @@ class RadicalTessellationContactConstruction const SimplePoint& ac_plane_center, const SimplePoint& ac_plane_normal, const UnsignedInt c_id, - Contour& contour) + Contour& contour) noexcept { UnsignedInt outsiders_count=0; for(Contour::iterator it=contour.begin();it!=contour.end();++it) @@ -347,7 +352,7 @@ class RadicalTessellationContactConstruction const SimplePoint& ac_plane_center, const SimplePoint& ac_plane_normal, const UnsignedInt c_id, - Contour& contour) + Contour& contour) noexcept { if(contour.size()<3) { @@ -431,7 +436,7 @@ class RadicalTessellationContactConstruction } } - static bool test_if_contour_is_still_cuttable(const SimplePoint& a_center, const SimplePoint& closest_possible_cut_point, const Contour& contour) + static bool test_if_contour_is_still_cuttable(const SimplePoint& a_center, const SimplePoint& closest_possible_cut_point, const Contour& contour) noexcept { bool cuttable=false; const Float dist_threshold=squared_distance_from_point_to_point(a_center, closest_possible_cut_point); @@ -447,7 +452,7 @@ class RadicalTessellationContactConstruction const SimplePoint& ic_axis, const UnsignedInt a_id, Contour& contour, - Float& sum_of_arc_angles) + Float& sum_of_arc_angles) noexcept { UnsignedInt outsiders_count=0; for(UnsignedInt i=0;i2 && equal(sum_of_arc_angles, PIVALUE*FLOATCONST(2.0), 0.001))) + { + sum_of_arc_angles=FLOATCONST(2.0); + contour.clear(); + } } } return (outsiders_count>0); } - static Float calculate_contour_area(const SimpleSphere& ic_sphere, const Contour& contour, SimplePoint& contour_barycenter) + static Float calculate_contour_area(const SimpleSphere& ic_sphere, const Contour& contour, SimplePoint& contour_barycenter) noexcept { Float area=FLOATCONST(0.0); @@ -605,7 +615,7 @@ class RadicalTessellationContactConstruction return area; } - static Float calculate_contour_solid_angle(const SimpleSphere& a, const SimpleSphere& b, const SimpleSphere& ic_sphere, const Contour& contour) + static Float calculate_contour_solid_angle(const SimpleSphere& a, const SimpleSphere& b, const SimpleSphere& ic_sphere, const Contour& contour) noexcept { Float turn_angle=FLOATCONST(0.0); @@ -657,16 +667,25 @@ class RadicalTessellationContactConstruction return solid_angle; } - static bool check_if_contour_is_central(const SimplePoint& center, const Contour& contour, const SimplePoint& contour_barycenter) + static bool check_if_contour_is_central(const SimplePoint& center, const Contour& contour, const SimplePoint& contour_barycenter) noexcept { - bool central=true; - for(UnsignedInt i=0;i -#include "preparation_for_tessellation.h" +#include "spheres_container.h" #include "simplified_aw_tessellation_contact_construction.h" #include "time_recorder.h" @@ -21,7 +21,7 @@ class SimplifiedAWTessellation UnsignedInt id_a; UnsignedInt id_b; - ContactDescriptorSummary() : + ContactDescriptorSummary() noexcept : area(FLOATCONST(0.0)), arc_length(FLOATCONST(0.0)), distance(FLOATCONST(0.0)), @@ -30,7 +30,7 @@ class SimplifiedAWTessellation { } - void set(const SimplifiedAWTessellationContactConstruction::ContactDescriptor& cd) + void set(const SimplifiedAWTessellationContactConstruction::ContactDescriptor& cd) noexcept { if(cd.area>FLOATCONST(0.0)) { @@ -42,7 +42,7 @@ class SimplifiedAWTessellation } } - void ensure_ids_ordered() + void ensure_ids_ordered() noexcept { if(id_a>id_b) { @@ -58,7 +58,7 @@ class SimplifiedAWTessellation Float distance; UnsignedInt count; - TotalContactDescriptorsSummary() : + TotalContactDescriptorsSummary() noexcept : area(FLOATCONST(0.0)), arc_length(FLOATCONST(0.0)), distance(FLOATCONST(-1.0)), @@ -66,7 +66,7 @@ class SimplifiedAWTessellation { } - void add(const ContactDescriptorSummary& cds) + void add(const ContactDescriptorSummary& cds) noexcept { if(cds.area>FLOATCONST(0.0)) { @@ -86,7 +86,7 @@ class SimplifiedAWTessellation std::vector contacts_summaries; TotalContactDescriptorsSummary total_contacts_summary; - Result() : total_spheres(0), total_collisions(0), total_relevant_collisions(0) + Result() noexcept : total_spheres(0), total_collisions(0), total_relevant_collisions(0) { } }; @@ -104,7 +104,7 @@ class SimplifiedAWTessellation static void construct_full_tessellation( const std::vector& spheres, - Result& result) + Result& result) noexcept { TimeRecorder time_recorder; ResultGraphics result_graphics; @@ -117,18 +117,21 @@ class SimplifiedAWTessellation const bool with_graphics, Result& result, ResultGraphics& result_graphics, - TimeRecorder& time_recorder) + TimeRecorder& time_recorder) noexcept { - PreparationForTessellation::Result preparation_result; - PreparationForTessellation::prepare_for_tessellation(spheres, grouping_of_spheres, preparation_result, time_recorder); + SpheresContainer spheres_container; + spheres_container.init(spheres, time_recorder); + + SpheresContainer::ResultOfPreparationForTessellation preparation_result; + spheres_container.prepare_for_tessellation(grouping_of_spheres, preparation_result, time_recorder); time_recorder.reset(); result=Result(); result_graphics=ResultGraphics(); - result.total_spheres=preparation_result.total_spheres; - result.total_collisions=preparation_result.total_collisions; + result.total_spheres=spheres_container.input_spheres().size(); + result.total_collisions=spheres_container.total_collisions(); result.total_relevant_collisions=preparation_result.relevant_collision_ids.size(); std::vector possible_contacts_summaries(preparation_result.relevant_collision_ids.size()); @@ -141,17 +144,21 @@ class SimplifiedAWTessellation time_recorder.record_elapsed_miliseconds_and_reset("allocate possible contact summaries"); - #pragma omp parallel +#ifdef VORONOTALT_OPENMP +#pragma omp parallel +#endif { SimplifiedAWTessellationContactConstruction::ContactDescriptor cd; cd.neighbor_descriptors.reserve(24); - #pragma omp for +#ifdef VORONOTALT_OPENMP +#pragma omp for +#endif for(UnsignedInt i=0;i& grouping_of_spheres, - GroupedResult& grouped_result) + GroupedResult& grouped_result) noexcept { TimeRecorder time_recorder; return group_results(full_result, grouping_of_spheres, grouped_result, time_recorder); @@ -235,7 +243,7 @@ class SimplifiedAWTessellation const Result& full_result, const std::vector& grouping_of_spheres, GroupedResult& grouped_result, - TimeRecorder& time_recorder) + TimeRecorder& time_recorder) noexcept { time_recorder.reset(); diff --git a/src/voronotalt/simplified_aw_tessellation_contact_construction.h b/src/voronotalt/simplified_aw_tessellation_contact_construction.h index 7c01fd79f..973661e1b 100644 --- a/src/voronotalt/simplified_aw_tessellation_contact_construction.h +++ b/src/voronotalt/simplified_aw_tessellation_contact_construction.h @@ -19,7 +19,7 @@ class SimplifiedAWTessellationContactConstruction UnsignedInt left_id; UnsignedInt right_id; - ContourPoint(const SimplePoint& p, const UnsignedInt left_id, const UnsignedInt right_id) : p(p), left_id(left_id), right_id(right_id) + ContourPoint(const SimplePoint& p, const UnsignedInt left_id, const UnsignedInt right_id) noexcept : p(p), left_id(left_id), right_id(right_id) { } }; @@ -31,11 +31,11 @@ class SimplifiedAWTessellationContactConstruction Float sort_value; UnsignedInt neighbor_id; - NeighborDescriptor() : sort_value(FLOATCONST(0.0)), neighbor_id(0) + NeighborDescriptor() noexcept : sort_value(FLOATCONST(0.0)), neighbor_id(0) { } - bool operator<(const NeighborDescriptor& d) const + bool operator<(const NeighborDescriptor& d) const noexcept { return (sort_value contours_graphics; - ContactDescriptorGraphics() + ContactDescriptorGraphics() noexcept { } - void clear() + void clear() noexcept { contours_graphics.clear(); } @@ -74,7 +74,7 @@ class SimplifiedAWTessellationContactConstruction UnsignedInt id_a; UnsignedInt id_b; - ContactDescriptor() : + ContactDescriptor() noexcept : area(FLOATCONST(0.0)), arc_length(FLOATCONST(0.0)), distance(FLOATCONST(0.0)), @@ -83,7 +83,7 @@ class SimplifiedAWTessellationContactConstruction { } - void clear() + void clear() noexcept { id_a=0; id_b=0; @@ -101,11 +101,11 @@ class SimplifiedAWTessellationContactConstruction const std::vector& spheres_exclusion_statuses, const UnsignedInt a_id, const UnsignedInt b_id, - const std::vector& a_neighbor_collisions, + const std::vector& a_neighbor_collisions, const Float step, const int projections, const bool record_graphics, - ContactDescriptor& result_contact_descriptor) + ContactDescriptor& result_contact_descriptor) noexcept { result_contact_descriptor.clear(); if(a_id=spheres_exclusion_statuses.size() || spheres_exclusion_statuses[neighbor_id]==0)) { const SimpleSphere& c=spheres[neighbor_id]; @@ -302,7 +302,7 @@ class SimplifiedAWTessellationContactConstruction class HyperboloidBetweenTwoSpheres { public: - static inline SimplePoint project_point_on_hyperboloid(const SimplePoint& p, const SimpleSphere& s1, const SimpleSphere& s2) + static inline SimplePoint project_point_on_hyperboloid(const SimplePoint& p, const SimpleSphere& s1, const SimpleSphere& s2) noexcept { if(s1.r>s2.r) { @@ -321,7 +321,7 @@ class SimplifiedAWTessellationContactConstruction } } - static inline Float intersect_vector_with_hyperboloid(const SimplePoint& a, const SimplePoint& b, const SimpleSphere& s1, const SimpleSphere& s2) + static inline Float intersect_vector_with_hyperboloid(const SimplePoint& a, const SimplePoint& b, const SimpleSphere& s1, const SimpleSphere& s2) noexcept { if(s1.r>s2.r) { @@ -350,7 +350,7 @@ class SimplifiedAWTessellationContactConstruction } private: - static inline Float project_point_on_simple_hyperboloid(const Float x, const Float y, const Float d, const Float r1, const Float r2) + static inline Float project_point_on_simple_hyperboloid(const Float x, const Float y, const Float d, const Float r1, const Float r2) noexcept { if(r1>r2) { @@ -363,7 +363,7 @@ class SimplifiedAWTessellationContactConstruction } } - static inline Float intersect_vector_with_simple_hyperboloid(const SimplePoint& a, const SimplePoint& b, const Float d, const Float r1, const Float r2) + static inline Float intersect_vector_with_simple_hyperboloid(const SimplePoint& a, const SimplePoint& b, const Float d, const Float r1, const Float r2) noexcept { if(r1>r2) { @@ -411,7 +411,7 @@ class SimplifiedAWTessellationContactConstruction const SimpleSphere& base, const SimplePoint& axis, const Float length_step, - Contour& result) + Contour& result) noexcept { const Float angle_step=std::max(std::min(length_step/base.r, PIVALUE/FLOATCONST(3.0)), PIVALUE/FLOATCONST(36.0)); const SimplePoint first_point=point_and_number_product(any_normal_of_vector(axis), base.r); @@ -427,7 +427,7 @@ class SimplifiedAWTessellationContactConstruction const SimpleSphere& c, const UnsignedInt c_id, Contour& contour, - std::list& segments) + std::list& segments) noexcept { const UnsignedInt outsiders_count=mark_contour(a, c, c_id, contour); if(outsiders_count>0) @@ -458,7 +458,7 @@ class SimplifiedAWTessellationContactConstruction const SimpleSphere& a, const SimpleSphere& c, const UnsignedInt c_id, - Contour& contour) + Contour& contour) noexcept { UnsignedInt outsiders_count=0; for(Contour::iterator it=contour.begin();it!=contour.end();++it) @@ -478,7 +478,7 @@ class SimplifiedAWTessellationContactConstruction const SimpleSphere& c, const UnsignedInt c_id, Contour& contour, - std::list& cuts) + std::list& cuts) noexcept { UnsignedInt cuts_count=0; Contour::iterator it=contour.begin(); @@ -517,7 +517,7 @@ class SimplifiedAWTessellationContactConstruction return cuts_count; } - static void order_cuts(std::list& cuts) + static void order_cuts(std::list& cuts) noexcept { Float sums[2]={FLOATCONST(0.0), FLOATCONST(0.0)}; for(int i=0;i<2;i++) @@ -552,7 +552,7 @@ class SimplifiedAWTessellationContactConstruction static UnsignedInt split_contour( Contour& contour, const std::list& ordered_cuts, - std::list& segments) + std::list& segments) noexcept { UnsignedInt segments_count=0; std::list::const_iterator it=ordered_cuts.begin(); @@ -597,7 +597,7 @@ class SimplifiedAWTessellationContactConstruction const UnsignedInt c_id, const Float step, const int projections, - Contour& contour) + Contour& contour) noexcept { Contour::iterator it=contour.begin(); while(it!=contour.end()) @@ -642,7 +642,7 @@ class SimplifiedAWTessellationContactConstruction } } - static bool check_contour_intersects_sphere(const SimpleSphere& shell, const Contour& contour) + static bool check_contour_intersects_sphere(const SimpleSphere& shell, const Contour& contour) noexcept { for(Contour::const_iterator it=contour.begin();it!=contour.end();++it) { @@ -654,7 +654,7 @@ class SimplifiedAWTessellationContactConstruction return false; } - static void filter_contours_intersecting_sphere(const SimpleSphere& shell, std::list& contours) + static void filter_contours_intersecting_sphere(const SimpleSphere& shell, std::list& contours) noexcept { std::list::iterator it=contours.begin(); while(it!=contours.end()) @@ -671,7 +671,7 @@ class SimplifiedAWTessellationContactConstruction } template - static Iterator get_left_iterator(List& container, const Iterator& iterator) + static Iterator get_left_iterator(List& container, const Iterator& iterator) noexcept { Iterator left_it=iterator; if(left_it==container.begin()) @@ -683,7 +683,7 @@ class SimplifiedAWTessellationContactConstruction } template - static Iterator get_right_iterator(List& container, const Iterator& iterator) + static Iterator get_right_iterator(List& container, const Iterator& iterator) noexcept { Iterator right_it=iterator; ++right_it; @@ -695,7 +695,7 @@ class SimplifiedAWTessellationContactConstruction } template - static void shift_list(List& list, const bool reverse) + static void shift_list(List& list, const bool reverse) noexcept { if(!reverse) { @@ -709,7 +709,7 @@ class SimplifiedAWTessellationContactConstruction } } - static bool collect_points_from_contour(const Contour& contour, std::vector& contour_points) + static bool collect_points_from_contour(const Contour& contour, std::vector& contour_points) noexcept { contour_points.clear(); if(!contour.empty()) @@ -723,7 +723,7 @@ class SimplifiedAWTessellationContactConstruction return (!contour_points.empty()); } - static Float calculate_area_from_contour(const Contour& contour, const SimpleSphere& sphere1, const SimpleSphere& sphere2, SimplePoint& contour_barycenter) + static Float calculate_area_from_contour(const Contour& contour, const SimpleSphere& sphere1, const SimpleSphere& sphere2, SimplePoint& contour_barycenter) noexcept { Float area=FLOATCONST(0.0); if(!contour.empty()) @@ -757,7 +757,7 @@ class SimplifiedAWTessellationContactConstruction return area; } - static Float calculate_arc_length_from_contour(const UnsignedInt a_id, const Contour& contour) + static Float calculate_arc_length_from_contour(const UnsignedInt a_id, const Contour& contour) noexcept { Float arc_length=FLOATCONST(0.0); for(Contour::const_iterator jt=contour.begin();jt!=contour.end();++jt) diff --git a/src/voronotalt/spheres_container.h b/src/voronotalt/spheres_container.h new file mode 100644 index 000000000..b01659e66 --- /dev/null +++ b/src/voronotalt/spheres_container.h @@ -0,0 +1,664 @@ +#ifndef VORONOTALT_SPHERES_CONTAINER_H_ +#define VORONOTALT_SPHERES_CONTAINER_H_ + +#include "spheres_searcher.h" +#include "periodic_box.h" +#include "time_recorder.h" + +namespace voronotalt +{ + +class SpheresContainer +{ +public: + struct ResultOfPreparationForTessellation + { + std::vector< std::pair > relevant_collision_ids; + + ResultOfPreparationForTessellation() noexcept + { + } + }; + + SpheresContainer() noexcept : total_collisions_(0) + { + } + + void init(const std::vector& input_spheres, TimeRecorder& time_recorder) noexcept + { + init(input_spheres, PeriodicBox(), time_recorder); + } + + void init(const std::vector& input_spheres, const PeriodicBox& periodic_box, TimeRecorder& time_recorder) noexcept + { + time_recorder.reset(); + + periodic_box_=periodic_box; + input_spheres_=input_spheres; + + if(periodic_box_.enabled()) + { + populated_spheres_.resize(input_spheres_.size()*27); + std::vector collected_indices; + for(UnsignedInt i=0;i& new_input_spheres, + const std::vector& provided_ids_of_changed_input_spheres, + const bool trust_provided_ids_of_changed_input_spheres, + std::vector& ids_of_changed_input_spheres, + std::vector& ids_of_affected_input_spheres, + TimeRecorder& time_recorder) noexcept + { + time_recorder.reset(); + + ids_of_changed_input_spheres.clear(); + ids_of_affected_input_spheres.clear(); + + if(new_input_spheres.size()!=input_spheres_.size()) + { + reinit(new_input_spheres, ids_of_changed_input_spheres, ids_of_affected_input_spheres, time_recorder); + return true; + } + + if(trust_provided_ids_of_changed_input_spheres) + { + ids_of_changed_input_spheres=provided_ids_of_changed_input_spheres; + } + else + { + for(UnsignedInt i=0;isize_threshold_for_full_reinit()) + { + reinit(new_input_spheres, ids_of_changed_input_spheres, ids_of_affected_input_spheres, time_recorder); + return true; + } + + for(UnsignedInt i=0;i=input_spheres_.size()) + { + reinit(new_input_spheres, ids_of_changed_input_spheres, ids_of_affected_input_spheres, time_recorder); + return true; + } + } + + { + bool update_needed=false; + for(UnsignedInt i=0;!update_needed && i::iterator it=std::lower_bound(ids_of_affected_input_spheres.begin(), ids_of_affected_input_spheres.end(), affected_sphere_id); + if(it==ids_of_affected_input_spheres.end() || (*it)!=affected_sphere_id) + { + ids_of_affected_input_spheres.insert(it, affected_sphere_id); + } + } + else + { + reinit(new_input_spheres, ids_of_changed_input_spheres, ids_of_affected_input_spheres, time_recorder); + return true; + } + } + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("gather affected spheres ids for update"); + + { + if(periodic_box_.enabled()) + { + std::vector ids_of_changed_populated_spheres; + ids_of_changed_populated_spheres.reserve(ids_of_changed_input_spheres.size()*27); + for(UnsignedInt i=0;i::iterator it=std::lower_bound(buffered_temporary_storage_.more_ids_of_affected_input_spheres.begin(), buffered_temporary_storage_.more_ids_of_affected_input_spheres.end(), affected_sphere_id); + if(it==buffered_temporary_storage_.more_ids_of_affected_input_spheres.end() || (*it)!=affected_sphere_id) + { + if(!std::binary_search(ids_of_affected_input_spheres.begin(), ids_of_affected_input_spheres.end(), affected_sphere_id)) + { + buffered_temporary_storage_.more_ids_of_affected_input_spheres.insert(it, affected_sphere_id); + } + } + } + } + + if(!buffered_temporary_storage_.more_ids_of_affected_input_spheres.empty()) + { + { +#ifdef VORONOTALT_OPENMP +#pragma omp parallel for +#endif + for(UnsignedInt i=0;i=input_spheres_.size() || id_of_masked_input_sphere>=all_exclusion_statuses_.size() || (new_exclusion_status ? all_exclusion_statuses_[id_of_masked_input_sphere]>0 : all_exclusion_statuses_[id_of_masked_input_sphere]<1)) + { + return false; + } + + all_exclusion_statuses_[id_of_masked_input_sphere]=(new_exclusion_status ? 1 : 0); + + if(periodic_box_.enabled()) + { + set_exclusion_status_periodic_instances(id_of_masked_input_sphere); + } + + return true; + } + + bool update_by_setting_exclusion_status(const UnsignedInt id_of_masked_input_sphere, const bool new_exclusion_status, std::vector& ids_of_changed_input_spheres, std::vector& ids_of_affected_input_spheres) noexcept + { + ids_of_changed_input_spheres.clear(); + ids_of_affected_input_spheres.clear(); + + if(id_of_masked_input_sphere>=input_spheres_.size() || id_of_masked_input_sphere>=all_exclusion_statuses_.size() || (new_exclusion_status ? all_exclusion_statuses_[id_of_masked_input_sphere]>0 : all_exclusion_statuses_[id_of_masked_input_sphere]<1)) + { + return false; + } + + ids_of_changed_input_spheres.push_back(id_of_masked_input_sphere); + ids_of_affected_input_spheres.push_back(id_of_masked_input_sphere); + + for(std::size_t j=0;j::iterator it=std::lower_bound(ids_of_affected_input_spheres.begin(), ids_of_affected_input_spheres.end(), affected_sphere_id); + if(it==ids_of_affected_input_spheres.end() || (*it)!=affected_sphere_id) + { + ids_of_affected_input_spheres.insert(it, affected_sphere_id); + } + } + + all_exclusion_statuses_[id_of_masked_input_sphere]=(new_exclusion_status ? 1 : 0); + + if(periodic_box_.enabled()) + { + set_exclusion_status_periodic_instances(id_of_masked_input_sphere); + } + + return true; + } + + void assign(const SpheresContainer& obj) noexcept + { + periodic_box_=obj.periodic_box_; + total_collisions_=obj.total_collisions_; + + spheres_searcher_.assign(obj.spheres_searcher_); + + input_spheres_.resize(obj.input_spheres_.size()); + populated_spheres_.resize(obj.populated_spheres_.size()); + all_exclusion_statuses_.resize(obj.all_exclusion_statuses_.size()); + all_colliding_ids_.resize(obj.all_colliding_ids_.size()); + + { +#ifdef VORONOTALT_OPENMP +#pragma omp parallel for +#endif + for(UnsignedInt i=0;i& subset_of_ids) noexcept + { + if(subset_of_ids.empty() + || obj.input_spheres_.empty() + || !periodic_box_.equals(obj.periodic_box_) + || input_spheres_.size()!=obj.input_spheres_.size() + || populated_spheres_.size()!=obj.populated_spheres_.size() + || all_exclusion_statuses_.size()!=obj.all_exclusion_statuses_.size() + || all_colliding_ids_.size()!=obj.all_colliding_ids_.size() + || subset_of_ids.size()>=size_threshold_for_full_reinit()) + { + assign(obj); + } + else + { + for(UnsignedInt i=0;i=input_spheres_.size()) + { + assign(obj); + return; + } + } + + periodic_box_=obj.periodic_box_; + total_collisions_=obj.total_collisions_; + + spheres_searcher_.assign(obj.spheres_searcher_); + + { +#ifdef VORONOTALT_OPENMP +#pragma omp parallel for +#endif + for(UnsignedInt i=0;i& input_spheres() const noexcept + { + return input_spheres_; + } + + const std::vector& populated_spheres() const noexcept + { + return populated_spheres_; + } + + const std::vector& all_exclusion_statuses() const noexcept + { + return all_exclusion_statuses_; + } + + const std::vector< std::vector >& all_colliding_ids() const noexcept + { + return all_colliding_ids_; + } + + UnsignedInt total_collisions() const noexcept + { + return total_collisions_; + } + + bool prepare_for_tessellation(const std::vector& grouping_of_spheres, ResultOfPreparationForTessellation& result, TimeRecorder& time_recorder) const noexcept + { + return prepare_for_tessellation(std::vector(), grouping_of_spheres, result, time_recorder); + } + + bool prepare_for_tessellation(const std::vector& involvement_of_spheres, const std::vector& grouping_of_spheres, ResultOfPreparationForTessellation& result, TimeRecorder& time_recorder) const noexcept + { + time_recorder.reset(); + + result.relevant_collision_ids.clear(); + + if(populated_spheres_.empty()) + { + return false; + } + + result.relevant_collision_ids.reserve(total_collisions_); + + for(UnsignedInt id_a=0;id_a=involvement_of_spheres.size() || involvement_of_spheres[id_a]>0) && all_exclusion_statuses_[id_a]==0) + { + for(UnsignedInt j=0;j=involvement_of_spheres.size() || involvement_of_spheres[id_b_canonical]>0) && all_exclusion_statuses_[id_b_canonical]==0) + { + if(id_b>=input_spheres_.size() || (all_colliding_ids_[id_a].size()=grouping_of_spheres.size() || id_b_canonical>=grouping_of_spheres.size() || grouping_of_spheres[id_a]!=grouping_of_spheres[id_b_canonical]) + { + result.relevant_collision_ids.push_back(std::pair(id_a, id_b)); + } + } + } + } + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("collect relevant collision indices"); + + return true; + } + +private: + struct BufferedTemporaryStorage + { + std::vector more_ids_of_affected_input_spheres; + std::vector merged_ids_of_affected_input_spheres; + + void clear() noexcept + { + more_ids_of_affected_input_spheres.clear(); + merged_ids_of_affected_input_spheres.clear(); + } + }; + + UnsignedInt size_threshold_for_full_reinit() const noexcept + { + return static_cast(input_spheres_.size()/2); + } + + void reinit(const std::vector& new_input_spheres, std::vector& ids_of_changed_input_spheres, std::vector& ids_of_affected_input_spheres, TimeRecorder& time_recorder) noexcept + { + init(new_input_spheres, periodic_box_, time_recorder); + ids_of_changed_input_spheres.clear(); + ids_of_affected_input_spheres.clear(); + } + + void set_sphere_periodic_instances(const UnsignedInt i, const bool collect_indices, std::vector& collected_indices) noexcept + { + if(i(sx), static_cast(sy), static_cast(sz)); + if(collect_indices) + { + collected_indices.push_back(mi); + } + g++; + } + } + } + } + } + } + } + + void set_exclusion_status_periodic_instances(const UnsignedInt i) noexcept + { + if(i input_spheres_; + std::vector populated_spheres_; + std::vector all_exclusion_statuses_; + SpheresSearcher spheres_searcher_; + std::vector< std::vector > all_colliding_ids_; + UnsignedInt total_collisions_; + BufferedTemporaryStorage buffered_temporary_storage_; +}; + +} + +#endif /* VORONOTALT_SPHERES_CONTAINER_H_ */ diff --git a/src/voronotalt/spheres_searcher.h b/src/voronotalt/spheres_searcher.h index 52a73a593..b34652854 100644 --- a/src/voronotalt/spheres_searcher.h +++ b/src/voronotalt/spheres_searcher.h @@ -2,6 +2,7 @@ #define VORONOTALT_SPHERES_SEARCHER_H_ #include +#include #include "basic_types_and_functions.h" @@ -11,71 +12,99 @@ namespace voronotalt class SpheresSearcher { public: - explicit SpheresSearcher(const std::vector& spheres) : spheres_(spheres), box_size_(FLOATCONST(1.0)) + SpheresSearcher() noexcept { - box_size_=calculate_grid_box_size(spheres_, box_size_); + } - for(UnsignedInt i=0;i& spheres) noexcept + { + init(spheres); + } + + void init(const std::vector& spheres) noexcept + { + spheres_=spheres; + grid_parameters_=GridParameters(spheres_); + init_boxes(); + } + + bool update(const std::vector& spheres, const std::vector& ids_to_update) noexcept + { + if(ids_to_update.empty()) { - const GridPoint gp(spheres_[i], box_size_); - if(i==0) - { - grid_offset_=gp; - grid_size_=gp; - } - else + return false; + } + + if(ids_to_update.size()>size_threshold_for_full_rebuild()) + { + init(spheres); + return true; + } + + for(UnsignedInt i=0;i=spheres.size() || !update_sphere(ids_to_update[i], spheres[ids_to_update[i]])) { - grid_offset_.x=std::min(grid_offset_.x, gp.x); - grid_offset_.y=std::min(grid_offset_.y, gp.y); - grid_offset_.z=std::min(grid_offset_.z, gp.z); - grid_size_.x=std::max(grid_size_.x, gp.x); - grid_size_.y=std::max(grid_size_.y, gp.y); - grid_size_.z=std::max(grid_size_.z, gp.z); + init(spheres); + return true; } } - grid_size_.x=grid_size_.x-grid_offset_.x+1; - grid_size_.y=grid_size_.y-grid_offset_.y+1; - grid_size_.z=grid_size_.z-grid_offset_.z+1; + return true; + } + + void assign(const SpheresSearcher& obj) noexcept + { + grid_parameters_=obj.grid_parameters_; - map_of_boxes_.resize(grid_size_.x*grid_size_.y*grid_size_.z, -1); + spheres_.resize(obj.spheres_.size()); + map_of_boxes_.resize(obj.map_of_boxes_.size()); + boxes_.resize(obj.boxes_.size()); - for(UnsignedInt i=0;i(boxes_.size()); - boxes_.push_back(std::vector(1, i)); + spheres_[i]=obj.spheres_[i]; } - else + } + + { +#ifdef VORONOTALT_OPENMP +#pragma omp parallel for +#endif + for(UnsignedInt i=0;i& spheres, const Float min_box_size) - { - Float box_size=min_box_size; - for(UnsignedInt i=0;i& colliding_ids, const bool discard_hidden, int& exclusion_status) const + const std::vector& searchable_spheres() const noexcept + { + return spheres_; + } + + bool find_colliding_ids(const UnsignedInt& central_id, std::vector& colliding_ids, const bool discard_hidden, int& exclusion_status) const noexcept { colliding_ids.clear(); exclusion_status=0; if(central_id=0) { const int box_id=map_of_boxes_[index]; @@ -110,7 +139,7 @@ class SpheresSearcher else if(!discard_hidden || !sphere_contains_sphere(central_sphere, candidate_sphere)) { - colliding_ids.push_back(id); + colliding_ids.push_back(ValuedID(distance_to_center_of_intersection_circle_of_two_spheres(central_sphere, candidate_sphere), id)); } } } @@ -120,6 +149,10 @@ class SpheresSearcher } } } + if(!colliding_ids.empty()) + { + std::sort(colliding_ids.begin(), colliding_ids.end()); + } return (!colliding_ids.empty()); } @@ -130,51 +163,205 @@ class SpheresSearcher int y; int z; - GridPoint() : x(0), y(0), z(0) + GridPoint() noexcept: x(0), y(0), z(0) { } - GridPoint(const SimpleSphere& s, const Float grid_step) + GridPoint(const SimpleSphere& s, const Float grid_step) noexcept { init(s, grid_step); } - GridPoint(const SimpleSphere& s, const Float grid_step, const GridPoint& grid_offset) + GridPoint(const SimpleSphere& s, const Float grid_step, const GridPoint& grid_offset) noexcept { init(s, grid_step, grid_offset); } - void init(const SimpleSphere& s, const Float grid_step) + void init(const SimpleSphere& s, const Float grid_step) noexcept { x=static_cast(s.p.x/grid_step); y=static_cast(s.p.y/grid_step); z=static_cast(s.p.z/grid_step); } - void init(const SimpleSphere& s, const Float grid_step, const GridPoint& grid_offset) + void init(const SimpleSphere& s, const Float grid_step, const GridPoint& grid_offset) noexcept { x=static_cast(s.p.x/grid_step)-grid_offset.x; y=static_cast(s.p.y/grid_step)-grid_offset.y; z=static_cast(s.p.z/grid_step)-grid_offset.z; } - int index(const GridPoint& grid_size) const + int index(const GridPoint& grid_size) const noexcept { - return ((x>=0 && y>=0 && z>=0 && x=0 && y>=0 && z>=0 && x& spheres_; - GridPoint grid_offset_; - GridPoint grid_size_; + struct GridParameters + { + GridPoint grid_offset; + GridPoint grid_size; + Float box_size; + int padding_; + + GridParameters() noexcept : box_size(FLOATCONST(1.0)), padding_(1) + { + grid_size.x=1; + grid_size.y=1; + grid_size.z=1; + } + + explicit GridParameters(const std::vector& spheres) noexcept : box_size(FLOATCONST(1.0)), padding_(1) + { + init(spheres); + } + + GridParameters(const std::vector& spheres, const int padding) noexcept : box_size(FLOATCONST(1.0)), padding_(padding) + { + init(spheres); + } + + void init(const std::vector& spheres) noexcept + { + box_size=FLOATCONST(1.0); + padding_=std::max(0, padding_); + + for(UnsignedInt i=0;i(spheres_.size()/2); + } + + void init_boxes() noexcept + { + map_of_boxes_.clear(); + boxes_.clear(); + + map_of_boxes_.resize(grid_parameters_.grid_size.x*grid_parameters_.grid_size.y*grid_parameters_.grid_size.z, -1); + + for(UnsignedInt i=0;i(boxes_.size()); + boxes_.push_back(std::vector(1, i)); + } + else + { + boxes_[box_id].push_back(i); + } + } + } + + bool update_sphere(const UnsignedInt& sphere_id, const SimpleSphere& moved_sphere) noexcept + { + if(sphere_id=0) + { + const GridPoint gp0(spheres_[sphere_id], grid_parameters_.box_size, grid_parameters_.grid_offset); + const int index0=gp0.index(grid_parameters_.grid_size); + if(index0>=0) + { + if(index1!=index0) + { + { + const int box_id0=map_of_boxes_[index0]; + if(box_id0<0) + { + return false; + } + else + { + std::vector& v0=boxes_[box_id0]; + std::vector::iterator it=std::lower_bound(v0.begin(), v0.end(), sphere_id); + if(it!=v0.end() && (*it)==sphere_id) + { + v0.erase(it); + } + } + } + + { + const int box_id1=map_of_boxes_[index1]; + if(box_id1<0) + { + map_of_boxes_[index1]=static_cast(boxes_.size()); + boxes_.push_back(std::vector(1, sphere_id)); + } + else + { + std::vector& v1=boxes_[box_id1]; + std::vector::iterator it=std::lower_bound(v1.begin(), v1.end(), sphere_id); + if(it==v1.end() || (*it)!=sphere_id) + { + v1.insert(it, sphere_id); + } + } + } + } + spheres_[sphere_id]=moved_sphere; + return true; + } + } + } + return false; + } + + std::vector spheres_; + GridParameters grid_parameters_; std::vector map_of_boxes_; std::vector< std::vector > boxes_; - Float box_size_; }; } diff --git a/src/voronotalt/time_recorder.h b/src/voronotalt/time_recorder.h index ac7be6cc7..683dc9d2b 100644 --- a/src/voronotalt/time_recorder.h +++ b/src/voronotalt/time_recorder.h @@ -11,15 +11,15 @@ class TimeRecorder { } - virtual ~TimeRecorder() + virtual ~TimeRecorder() noexcept { } - virtual void reset() + virtual void reset() noexcept { } - virtual void record_elapsed_miliseconds_and_reset(const char* /*message*/) + virtual void record_elapsed_miliseconds_and_reset(const char* /*message*/) noexcept { } }; diff --git a/src/voronotalt/updateable_radical_tessellation.h b/src/voronotalt/updateable_radical_tessellation.h new file mode 100644 index 000000000..853f81fcf --- /dev/null +++ b/src/voronotalt/updateable_radical_tessellation.h @@ -0,0 +1,671 @@ +#ifndef VORONOTALT_UPDATEABLE_RADICAL_TESSELLATION_H_ +#define VORONOTALT_UPDATEABLE_RADICAL_TESSELLATION_H_ + +#include "radical_tessellation.h" + +namespace voronotalt +{ + +class UpdateableRadicalTessellation +{ +public: + struct Result + { + std::vector cells_summaries; + std::vector< std::vector > contacts_summaries; + std::vector< std::vector > contacts_summaries_with_redundancy_in_periodic_box; + + Result() noexcept + { + } + + bool empty() noexcept + { + return (cells_summaries.empty() || contacts_summaries.empty()); + } + }; + + struct ResultSummary + { + RadicalTessellation::TotalContactDescriptorsSummary total_contacts_summary; + RadicalTessellation::TotalCellContactDescriptorsSummary total_cells_summary; + + ResultSummary() noexcept + { + } + }; + + UpdateableRadicalTessellation() noexcept : backup_enabled_(false), in_sync_with_backup_(false) + { + } + + explicit UpdateableRadicalTessellation(const bool backup_enabled) noexcept : backup_enabled_(backup_enabled), in_sync_with_backup_(false) + { + } + + bool init(const std::vector& input_spheres) noexcept + { + TimeRecorder time_recorder; + return init(input_spheres, time_recorder); + } + + bool init(const std::vector& input_spheres, TimeRecorder& time_recorder) noexcept + { + return init(input_spheres, PeriodicBox(), time_recorder); + } + + bool init(const std::vector& input_spheres, const PeriodicBox& periodic_box) noexcept + { + TimeRecorder time_recorder; + return init(input_spheres, periodic_box, time_recorder); + } + + bool init(const std::vector& input_spheres, const PeriodicBox& periodic_box, TimeRecorder& time_recorder) noexcept + { + prepare_for_possible_init_or_update(time_recorder); + + state_.spheres_container.init(input_spheres, periodic_box, time_recorder); + + in_sync_with_backup_=false; + + RadicalTessellation::ResultGraphics result_graphics; + RadicalTessellation::construct_full_tessellation(state_.spheres_container, std::vector(), std::vector(), false, true, FLOATCONST(0.0), std::vector(), buffered_temporary_storage_.tessellation_result, result_graphics, time_recorder); + + involvement_of_spheres_for_update_.clear(); + + init_result_from_tessellation_result(); + + return !state_.result.empty(); + } + + bool update(const std::vector& new_input_spheres) noexcept + { + TimeRecorder time_recorder; + return update(new_input_spheres, std::vector(), false, time_recorder); + } + + bool update(const std::vector& new_input_spheres, TimeRecorder& time_recorder) noexcept + { + return update(new_input_spheres, std::vector(), false, time_recorder); + } + + bool update(const std::vector& new_input_spheres, const std::vector& ids_of_changed_input_spheres) noexcept + { + TimeRecorder time_recorder; + return update(new_input_spheres, ids_of_changed_input_spheres, true, time_recorder); + } + + bool update(const std::vector& new_input_spheres, const std::vector& ids_of_changed_input_spheres, TimeRecorder& time_recorder) noexcept + { + return update(new_input_spheres, ids_of_changed_input_spheres, true, time_recorder); + } + + bool update(const std::vector& new_input_spheres, const std::vector& provided_ids_of_changed_input_spheres, const bool trust_provided_ids_of_changed_input_spheres, TimeRecorder& time_recorder) noexcept + { + prepare_for_possible_init_or_update(time_recorder); + + if(trust_provided_ids_of_changed_input_spheres && provided_ids_of_changed_input_spheres.empty()) + { + return false; + } + + if(!state_.spheres_container.update(new_input_spheres, provided_ids_of_changed_input_spheres, trust_provided_ids_of_changed_input_spheres, state_.ids_of_changed_input_spheres, state_.ids_of_affected_input_spheres, time_recorder)) + { + return false; + } + + in_sync_with_backup_=false; + + if(state_.ids_of_affected_input_spheres.empty()) + { + RadicalTessellation::ResultGraphics result_graphics; + RadicalTessellation::construct_full_tessellation(state_.spheres_container, std::vector(), std::vector(), false, true, FLOATCONST(0.0), std::vector(), buffered_temporary_storage_.tessellation_result, result_graphics, time_recorder); + init_result_from_tessellation_result(); + } + else + { + update_using_current_state(time_recorder); + } + + return true; + } + + bool update_by_setting_exclusion_mask(const UnsignedInt id_of_targeted_input_sphere, const bool new_exclusion_status) noexcept + { + TimeRecorder time_recorder; + return update_by_setting_exclusion_mask(id_of_targeted_input_sphere, new_exclusion_status, time_recorder); + } + + bool update_by_setting_exclusion_mask(const UnsignedInt id_of_targeted_input_sphere, const bool new_exclusion_status, TimeRecorder& time_recorder) noexcept + { + if(state_.result.empty() || id_of_targeted_input_sphere>=state_.result.contacts_summaries.size() + || id_of_targeted_input_sphere>=state_.spheres_container.all_exclusion_statuses().size() + || (new_exclusion_status ? state_.spheres_container.all_exclusion_statuses()[id_of_targeted_input_sphere]>0 : state_.spheres_container.all_exclusion_statuses()[id_of_targeted_input_sphere]<1)) + { + return false; + } + + prepare_for_possible_init_or_update(time_recorder); + + if(new_exclusion_status) + { + if(!state_.spheres_container.update_by_setting_exclusion_status(id_of_targeted_input_sphere, true)) + { + return false; + } + + state_.ids_of_changed_input_spheres.push_back(id_of_targeted_input_sphere); + state_.ids_of_affected_input_spheres.push_back(id_of_targeted_input_sphere); + + for(UnsignedInt i=0;i >& all_result_volumes_for_contacts_summaries) noexcept + { + all_result_volumes_for_contacts_summaries.clear(); + + if(!backup_enabled_ || state_.result.empty()) + { + return false; + } + + const UnsignedInt N=state_.result.contacts_summaries.size(); + + std::vector< std::vector > all_neighbor_ids(N); + std::vector< std::vector > all_volumes_after_masking(N); + + for(UnsignedInt i=0;i& i_neighbor_ids=all_neighbor_ids[i]; + + { + const std::vector& cdss=state_.result.contacts_summaries[i]; + i_neighbor_ids.resize(cdss.size()); + for(UnsignedInt j=0;j& i_volumes_after_masking=all_volumes_after_masking[i]; + i_volumes_after_masking.resize(i_neighbor_ids.size(), FLOATCONST(0.0)); + + if(update_by_setting_exclusion_mask(i, true)) + { + for(UnsignedInt j=0;j& input_spheres() const noexcept + { + return state_.spheres_container.input_spheres(); + } + + bool exclusion_status_of_input_sphere(const UnsignedInt id_of_input_sphere) const noexcept + { + return (id_of_input_sphere0); + } + + const Result& result() const noexcept + { + return state_.result; + } + + ResultSummary result_summary() const noexcept + { + ResultSummary result_summary; + for(UnsignedInt i=0;i& last_update_ids_of_changed_input_spheres() const noexcept + { + return state_.ids_of_changed_input_spheres; + } + + const std::vector& last_update_ids_of_affected_input_spheres() const noexcept + { + return state_.ids_of_affected_input_spheres; + } + + bool last_update_was_full_reinit() const noexcept + { + return state_.last_update_was_full_reinit; + } + +private: + struct BufferedTemporaryStorage + { + RadicalTessellation::Result tessellation_result; + + void clear() noexcept + { + tessellation_result.clear(); + } + }; + + class ConditionToRemoveContact + { + public: + explicit ConditionToRemoveContact(const std::vector& involvement) noexcept : involvement_(involvement) + { + } + + bool operator()(const RadicalTessellation::ContactDescriptorSummary& cds) noexcept + { + return (involvement_.empty() || (involvement_[cds.id_a%involvement_.size()]>0 && involvement_[cds.id_b%involvement_.size()]>0)); + } + + private: + const std::vector& involvement_; + }; + + struct State + { + SpheresContainer spheres_container; + Result result; + std::vector ids_of_changed_input_spheres; + std::vector ids_of_affected_input_spheres; + bool last_update_was_full_reinit; + + State() noexcept : last_update_was_full_reinit(true) + { + } + + void assign(const State& obj) noexcept + { + spheres_container.assign(obj.spheres_container); + + result.cells_summaries.resize(obj.result.cells_summaries.size()); + result.contacts_summaries.resize(obj.result.contacts_summaries.size()); + result.contacts_summaries_with_redundancy_in_periodic_box.resize(obj.result.contacts_summaries_with_redundancy_in_periodic_box.size()); + ids_of_changed_input_spheres.resize(obj.ids_of_changed_input_spheres.size()); + ids_of_affected_input_spheres.resize(obj.ids_of_affected_input_spheres.size()); + + { + #pragma omp parallel for + for(UnsignedInt i=0;i& subset_of_ids_of_spheres) noexcept + { + if(!assign_everything && subset_of_ids_of_spheres.empty()) + { + return; + } + + if(assign_everything + || result.cells_summaries.size()!=obj.result.cells_summaries.size() + || result.contacts_summaries.size()!=obj.result.contacts_summaries.size() + || result.contacts_summaries_with_redundancy_in_periodic_box.size()!=obj.result.contacts_summaries_with_redundancy_in_periodic_box.size()) + { + assign(obj); + return; + } + + const bool periodic=!obj.result.contacts_summaries_with_redundancy_in_periodic_box.empty(); + + for(UnsignedInt i=0;i=result.cells_summaries.size() || (periodic && sphere_id>=result.contacts_summaries_with_redundancy_in_periodic_box.size())) + { + assign(obj); + return; + } + } + + spheres_container.assign(obj.spheres_container, ids_of_changed_input_spheres); + + { + #pragma omp parallel for + for(UnsignedInt i=0;i& all_contacts_summaries=(buffered_temporary_storage_.tessellation_result.contacts_summaries_with_redundancy_in_periodic_box.empty() ? buffered_temporary_storage_.tessellation_result.contacts_summaries : buffered_temporary_storage_.tessellation_result.contacts_summaries_with_redundancy_in_periodic_box); + for(UnsignedInt i=0;i(), false, false, FLOATCONST(0.0), std::vector(), buffered_temporary_storage_.tessellation_result, result_graphics, time_recorder); + + { + const ConditionToRemoveContact condition_to_remove_contact(involvement_of_spheres_for_update_); + + for(UnsignedInt i=0;i& v=state_.result.contacts_summaries[sphere_id]; + std::vector::iterator it=std::remove_if(v.begin(), v.end(), condition_to_remove_contact); + v.erase(it, v.end()); + } + if(!state_.result.contacts_summaries_with_redundancy_in_periodic_box.empty()) + { + std::vector& v=state_.result.contacts_summaries_with_redundancy_in_periodic_box[sphere_id]; + std::vector::iterator it=std::remove_if(v.begin(), v.end(), condition_to_remove_contact); + v.erase(it, v.end()); + } + } + } + + if(!buffered_temporary_storage_.tessellation_result.contacts_summaries.empty()) + { + for(UnsignedInt i=0;i& all_contacts_summaries=(buffered_temporary_storage_.tessellation_result.contacts_summaries_with_redundancy_in_periodic_box.empty() ? buffered_temporary_storage_.tessellation_result.contacts_summaries : buffered_temporary_storage_.tessellation_result.contacts_summaries_with_redundancy_in_periodic_box); + + for(UnsignedInt i=0;i >& all_contacts_summaries=(state_.result.contacts_summaries_with_redundancy_in_periodic_box.empty() ? state_.result.contacts_summaries : state_.result.contacts_summaries_with_redundancy_in_periodic_box); + + for(UnsignedInt i=0;i& v=all_contacts_summaries[sphere_id]; + for(UnsignedInt j=0;j involvement_of_spheres_for_update_; + BufferedTemporaryStorage buffered_temporary_storage_; +}; + +} + +#endif /* VORONOTALT_UPDATEABLE_RADICAL_TESSELLATION_H_ */ diff --git a/src/voronotalt/voronotalt.h b/src/voronotalt/voronotalt.h index d8d651458..8ed18c8e7 100644 --- a/src/voronotalt/voronotalt.h +++ b/src/voronotalt/voronotalt.h @@ -4,5 +4,6 @@ #include "conversion_to_input.h" #include "radical_tessellation.h" #include "simplified_aw_tessellation.h" +#include "updateable_radical_tessellation.h" #endif /* VORONOTALT_VORONOTALT_H_ */