Merge pull request #2 from Ladme/v0.4

v0.4

Author: Ladme

CHANGELOG.md

@@ -1,3 +1,9 @@
+## Version 0.4.0
+- **Geometry selection:** Added the ability to select a geometric region for analysis. Users can now specify cuboidal, spherical, or cylindrical regions, and order parameters will be calculated only for bonds located within the selected region.
+- **Support for reading GRO, PDB, and PQR files:** These file formats are now supported as input structure files. In some cases, an additional "bonds" file specifying the system's connectivity may be required. Refer to the manual for more details.
+- **Manual assignment of lipids to leaflets:** Lipids can now be manually assigned to leaflets using a provided leaflet assignment file. Refer to the manual for detailed instructions.
+- **Calculating average results for the entire system:** YAML and TAB files now include information about the average order parameters calculated across all bonds and molecule types in the system. Additionally, ordermaps are generated for the entire system.
+
 ## Version 0.3.0
 - **Error estimation and convergence analysis**: Implemented error estimation and convergence analysis. Refer to the corresponding section of the manual for more details.
 - **Leaflet classification**: Leaflet classification can now also be performed either every N analyzed trajectory frames or only once at the start of the analysis.

Cargo.lock

@@ -1,8 +1,9 @@
 [package]
 name = "gorder"
-version = "0.3.0"
+version = "0.4.0"
 description = "Calculating lipid order parameters from Gromacs simulations"
 edition = "2021"
+rust-version = "1.82"
 license = "MIT"
 repository = "https://github.com/Ladme/gorder"
 keywords = ["gromacs", "molecular-dynamics", "order", "nmr", "lipids"]
@@ -11,7 +12,11 @@ exclude = ["/tests", "/validation"]
 
 [dependencies]
 backitup = "0.1.1"
-groan_rs = { version = "0.10.0-dev.1", features = ["parallel", "no-xdrfile"] }
+groan_rs = { version = "0.10.0-dev.2", features = [
+    "parallel",
+    "no-xdrfile",
+    "serde",
+] }
 colored = "2.1"
 thiserror = "2.0.3"
 derive_builder = "0.20.2"
@@ -32,3 +37,4 @@ parking_lot = "0.12.3"
 approx = "0.5.1"
 tempfile = "3.14.0"
 assert_cmd = "2.0.16"
+rand = "0.8.5"

Cargo.toml

@@ -20,8 +20,8 @@ $ cargo install gorder
 structure: system.tpr
 trajectory: md.xtc     # use your MD trajectory directly - no PBC handling or molecule fixing needed
 analysis_type: !AAOrder
-    heavy_atoms: "@membrane and element name carbon"
-    hydrogens: "@membrane and element name hydrogen"
+    heavy_atoms: "element name carbon"
+    hydrogens: "element name hydrogen"
 output: order.yaml
 ```
 
@@ -38,19 +38,19 @@ $ gorder YOUR_INPUT_YAML_FILE
 ## Features
 - **Atomistic and coarse-grained systems.** `gorder` is able to calculate atomistic and coarse-grained order parameters for individual bonds of individual lipid types.
 - **Powerful selection language.** `gorder` allows for simple yet powerful atom selection using a VMD-like selection language, supporting regular expressions and groups from NDX files. (Read more about the language [here](https://ladme.github.io/gsl-guide/).)
-- **Automatic identification of molecule types.** `gorder` automatically recognizes bonds and classifies molecule types based on their topology.
+- **Automatic identification of molecule types.** `gorder` automatically recognizes bonds and classifies molecule types based on their topology. Order parameters are calculated and reported separately for each molecule type.
 - **Various output formats.** `gorder` can output results in YAML, XVG, CSV, and custom "table" format.
+- **Supports any force-field.** `gorder` is completely force-field agnostic. Martini? CHARMM? Slipids? Your own toy force-field? As long as your lipids have bonds, it will work.
 - **Leaflet-wise analysis.** `gorder` can perform scrambling-safe assignment of lipids to membrane leaflets using three different methods, and then calculate lipid order parameters for individual leaflets.
 - **Order parameter maps.** `gorder` can construct 2D maps of order parameters, so you know what parts of the membrane are ordered and disordered.
 - **Error estimation.** `gorder` can automatically estimate the error of the analysis and indicate how well your analysis has converged.
-- **Supports any force-field.** `gorder` is completely force-field agnostic. Martini? CHARMM? Slipids? Your own toy force-field? As long as your lipids have bonds, it will work.
+- **Analysis of specific membrane regions.** `gorder` can dynamically select lipids in a specified part of the membrane and calculate order parameters only for them.
 - **Extremely fast.** `gorder` is extremely fast (see [below](#benchmarking)) due to its ability to read only the necessary atoms from XTC files and its support for multithreading.
 
 ## Planned
-- [ ] Dynamic selection of lipids for order parameter calculation based on geometric conditions (i.e., only calculating order parameters from a part of a membrane).
-- [ ] Dynamic membrane normal calculation, supporting membrane vesicles.
-- [ ] Python API: using `gorder` as a Python library.
 - [ ] United-atom order parameters.
+- [ ] Python API: using `gorder` as a Python library.
+- [ ] Dynamic membrane normal calculation, supporting membrane vesicles.
 - [ ] Improved multithreading (currently, multithreading is only implemented at the trajectory reading level).
 
 ## Validation

README.md

@@ -4,21 +4,24 @@
 //! Contains the implementation of the calculation of the atomistic order parameters.
 
 use super::{common::macros::group_name, topology::SystemTopology};
-use crate::analysis::common::{analyze_frame, prepare_master_group, sanity_check_molecules};
+use crate::analysis::common::{
+    analyze_frame, prepare_geometry_selection, prepare_master_group, sanity_check_molecules,
+};
+use crate::analysis::structure;
 use crate::errors::{AnalysisError, TopologyError};
+use crate::input::LeafletClassification;
 use crate::presentation::aaresults::AAOrderResults;
 use crate::presentation::{AnalysisResults, OrderResults};
 use crate::{input::Analysis, PANIC_MESSAGE};
 
+use groan_rs::prelude::ProgressPrinter;
 use groan_rs::prelude::{GroupXtcReader, OrderedAtomIterator};
-use groan_rs::{files::FileType, prelude::ProgressPrinter, system::System};
 
 /// Calculate the atomistic order parameters.
-pub(super) fn analyze_atomistic<'a>(
+pub(super) fn analyze_atomistic(
     analysis: Analysis,
 ) -> Result<AnalysisResults, Box<dyn std::error::Error + Send + Sync>> {
-    let mut system = System::from_file_with_format(analysis.structure(), FileType::TPR)?;
-    log::info!("Read molecular topology from '{}'.", analysis.structure());
+    let mut system = structure::read_structure_and_topology(&analysis)?;
 
     if let Some(ndx) = analysis.index() {
         system.read_ndx(ndx)?;
@@ -83,6 +86,9 @@ pub(super) fn analyze_atomistic<'a>(
         leaflet.prepare_system(&mut system)?;
     }
 
+    let geom = prepare_geometry_selection(analysis.geometry().as_ref(), &mut system)?;
+    geom.info();
+
     log::info!("Detecting molecule types...");
     log::logger().flush();
 
@@ -103,15 +109,22 @@ pub(super) fn analyze_atomistic<'a>(
         return Ok(AnalysisResults::AA(AAOrderResults::empty(analysis)));
     }
 
-    let data = SystemTopology::new(
+    let mut data = SystemTopology::new(
         molecules,
         analysis.membrane_normal().into(),
         analysis.estimate_error().clone(),
         analysis.step(),
         analysis.n_threads(),
+        geom,
     );
+
     data.info();
 
+    // finalize the manual leaflet classification
+    if let Some(LeafletClassification::Manual(params)) = analysis.leaflets() {
+        data.finalize_manual_leaflet_classification(params)?;
+    }
+
     let progress_printer = if analysis.silent() {
         None
     } else {
@@ -150,6 +163,9 @@ pub(super) fn analyze_atomistic<'a>(
         progress_printer,
     )?;
 
+    if let Some(LeafletClassification::Manual(_)) = analysis.leaflets() {
+        result.validate_manual_leaflet_classification(&analysis)?;
+    }
     result.log_total_analyzed_frames();
 
     // print basic info about error estimation
@@ -163,11 +179,14 @@ pub(super) fn analyze_atomistic<'a>(
 #[cfg(test)]
 mod tests {
     use approx::assert_relative_eq;
-    use groan_rs::prelude::Dimension;
+    use groan_rs::{prelude::Dimension, system::System};
 
     use super::*;
     use crate::{
-        analysis::molecule::{BondType, MoleculeType},
+        analysis::{
+            geometry::GeometrySelectionType,
+            molecule::{BondType, MoleculeType},
+        },
         input::leaflets::LeafletClassification,
     };
 
@@ -209,7 +228,14 @@ mod tests {
 
         (
             system,
-            SystemTopology::new(molecules, Dimension::Z, None, 1, 1),
+            SystemTopology::new(
+                molecules,
+                Dimension::Z,
+                None,
+                1,
+                1,
+                GeometrySelectionType::default(),
+            ),
         )
     }