groan_rs: 基于Rust的Gromacs分析库

Rust库，用于分析Gromacs模拟。目前处于开发早期阶段：任何内容都可能随时中断、更改或停止工作。

使用方法

运行

$ cargo add groan_rs

在Rust代码中导入该crate

use groan_rs::prelude::*;

许可证

此库根据MIT许可证发布。

完整文档

有关完整文档，请参阅相应的docs.rs页面。

lib.rs:

groan_rs: 基于Rust的Gromacs分析库

Rust库，用于分析Gromacs模拟。目前处于开发早期阶段：任何内容都可能随时中断、更改或停止工作。

它能做什么

• 读取和写入gro、pdb、pqr、ndx、xtc和trr文件。
• 从tpr文件中读取系统的拓扑和结构（使用minitprcrate）。
• 遍历原子并访问它们的属性，包括连接性（键）。
• 使用类似VMD的选择语言选择原子。
• 计算原子间的距离，同时遵守周期性边界条件。
• 根据几何条件选择原子。
• 分配元素到原子并猜测连接性。
• 计算任何原子组的几何中心或质量。
• 在一个模拟盒子中定位原子组。
• 通过提供实用数据结构（例如，GridMaps）帮助进行特定分析。
• 一些其他小功能...

它目前不能做什么

• 与非正交周期性边界条件一起工作。
• 无需额外工具即可进行结构和动力学的高级分析。（但是groan_rs库试图使实现自己的功能变得简单！)

使用方法

运行

$ cargo add groan_rs

在Rust代码中导入该crate

use groan_rs::prelude::*;

示例

分析结构文件

读取gro文件和ndx文件，并计算蛋白质的几何中心。

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file
    let mut system = System::from_file("system.gro")?;
    // `groan_rs` also supports pdb, pqr, and tpr files which
    // can be read using the same function as above

    // read an ndx file
    system.read_ndx("index.ndx")?;

    // calculate the center of geometry of a protein
    // note that the ndx file `index.ndx` must contain a group `Protein` for this to work
    let center = system.group_get_center("Protein")?;

    // print the result
    println!("{:?}", center);

    Ok(())
}

选择原子并创建组

读取gro文件，并选择属于系统特定部分的原子。groan_rs使用groan选择语言（GSL）来选择原子。GSL类似于VMD查询语言，见下文。

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file
    let mut system = System::from_file("system.gro")?;

    // select atoms using groan selection language creating a group 'My Group'
    system.group_create("My Group", "serial 25-28 or (resname POPC and name P)")?;
    // the group 'My Group' is now tightly associated with the system

    // we can now use the previously created group to construct another group
    // note that
    // a) each atom can be in any number of groups
    // b) the group names are case-sensitive
    system.group_create("my group", "'My Group' || resid 87 to 124")?;

    // we can then perform operations with the groups, e.g. write them into separate pdb files
    system.group_write_pdb("My Group", "My_Group.pdb", false)?;
    system.group_write_pdb("my group", "my_group.pdb", false)?;

    // each system also by default contains two groups consisting of all atoms in the system
    // these groups are called 'All' and 'all'
    assert!(system.group_exists("all"));
    assert!(system.group_exists("All"));

    // if you read an ndx file into the system like this:
    system.read_ndx("index.ndx")?;
    // the groups defined in the ndx file will be associated
    // with the system in the same way as the manually created groups

    Ok(())
}

分析轨迹文件

读取xtc文件，并计算从100纳秒开始到300纳秒结束的每一帧中两个原子组之间的距离。（groan_rs支持过程式和函数式方法。）

过程式方法

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file
    let mut system = System::from_file("structure.gro")?;

    // create the groups we are interested in
    // `groan_rs` uses VMD-like selection language for specifying groups of atoms
    system.group_create("group 1", "serial 1 to 5")?;
    system.group_create("group 2", "resid 45")?;

    // prepare a vector for the calculated distances
    let mut distances = Vec::new();

    // read the trajectory calculating the distance between the groups for each frame
    for frame in system
        .xtc_iter("files/md_short.xtc")?
        // we are only interested in frames between 100 and 300 ns
        .with_range(100_000.0, 300_000.)?
   {
        // check that the xtc frame has been read correctly
        let frame = frame?;
        // calculate the distance and put it into the vector
        distances.push(
            frame
            .group_distance("group 1", "group 2", Dimension::XYZ)
            .expect("Groups do not exist but they should."),
        );
    }

    // print the calculated distances
    println!("{:?}", distances);

    Ok(())
}

函数式方法

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file
    let mut system = System::from_file("structure.gro")?;

    // create the groups we are interested in
    // `groan_rs` uses VMD-like selection language for specifying groups of atoms
    system.group_create("group 1", "serial 1 to 5")?;
    system.group_create("group 2", "resid 45")?;

    // read the trajectory calculating distance between the groups
    // for each frame in the time range 100-300 ns
    // and collect the results into a vector
    let distances: Vec<f32> = system
        // read an xtc trajectory
        .xtc_iter("trajectory.xtc")?
        // we are only interested in frames between 100 and 300 ns
        .with_range(100_000.0, 300_000.)?
        // calculate distance between the groups for each frame
        .map(|frame| {
            // check that the xtc frame has been read correctly
            let frame = frame?;
            // calculate the distance
            Ok(frame
                .group_distance("group 1", "group 2", Dimension::XYZ)
                .expect("Groups do not exist but they should."))
        })
        // collect the calculated distances
        // if any error occured while reading the trajectory, propagate it
        .collect::<Result<Vec<f32>, Box<dyn Error + Send + Sync>>>()?;

    // print the calculated distances
    println!("{:?}", distances);

    Ok(())
}

请注意，with_range是一个非常有效的方法，它将跳过不在指定范围内的xtc帧，而无需实际读取这些帧中原子的属性。

您还可以让轨迹迭代器将轨迹读取信息打印到标准输出

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    let mut system = System::from_file("structure.gro")?;

    system.group_create("group 1", "serial 1 to 5")?;
    system.group_create("group 2", "resid 45")?;

    // create default progress printer
    let printer = ProgressPrinter::new();

    let distances: Vec<f32> = system
        .xtc_iter("trajectory.xtc")?
        // attach progress printer to the iterator
        .print_progress(printer)
        .map(|frame| {
            let frame = frame?;
            Ok(frame
                .group_distance("group 1", "group 2", Dimension::XYZ)
                .expect("Groups do not exist but they should."))
        })
        .collect::<Result<Vec<f32>, Box<dyn Error + Send + Sync>>>()?;

    println!("{:?}", distances);

    Ok(())
}

在trr和xtc文件之间转换

读取trr文件，并将粒子的位置写入新的xtc文件。

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file
    let mut system = System::from_file("structure.gro")?;

    // create an output xtc file
    // the `XtcWriter` structure is tightly coupled with the corresponding `System` structure
    // if the `System` is updated, the `XtcWriter` reflects this change
    let mut writer = XtcWriter::new(&system, "output.xtc")?;

    // iterate through the trr trajectory
    // the `trr_iter` (as well as `xtc_iter`) function updates the `System` structure
    // with which it is associated with the properties of the system in the current frame
    for frame in system.trr_iter("input.trr")? {
        // check that the trr frame has been read correctly
        frame?;
        // write the current frame into the output xtc file
        writer.write_frame()?;
    }

    Ok(())
}

遍历原子

遍历两个不同组中的原子，收集它们的信息或修改它们。

过程式方法

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file and an ndx file
    let mut system = System::from_file("system.gro")?;
    system.read_ndx("index.ndx")?;

    // immutably iterate over atoms of the group 'Protein' collecting atom names
    let mut names = Vec::new();
    for atom in system.group_iter("Protein")? {
        names.push(atom.get_atom_name().to_owned());
    }
    // (you can also iterate over all the atoms in the system using `System::atoms_iter`)

    // print the names
    println!("{:?}", names);

    // mutably iterate over atoms of the group 'Membrane' changing their residue names
    for atom in system.group_iter_mut("Membrane")? {
        atom.set_residue_name("MEMB");
    }
    // (you can also iterate over all the atoms in the system using `System::atoms_iter_mut`)

    Ok(())
}

函数式方法

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file and an ndx file
    let mut system = System::from_file("system.gro")?;
    system.read_ndx("index.ndx")?;

    // immutably iterate over atoms of the group 'Protein' collecting atom names
    let names: Vec<String> = system
        .group_iter("Protein")?
        .map(|atom| atom.get_atom_name().to_owned())
        .collect();
    // (you can also iterate over all the atoms in the system using `System::atoms_iter`)

    // print the names
    println!("{:?}", names);

    // mutably iterate over atoms of the group 'Membrane' changing their residue names
    system
        .group_iter_mut("Membrane")?
        .for_each(|atom| atom.set_residue_name("MEMB"));
    // (you can also iterate over all the atoms in the system using `System::atoms_iter_mut`)

    Ok(())
}

根据几何标准过滤和提取原子

过滤位于指定位置和尺寸的圆柱体内的原子，并将它们提取到单独的向量中。

过程式方法

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file
    let system = System::from_file("system.gro")?;

    // extract atoms that are located inside a specified cylinder
    let mut inside_cylinder = Vec::new();
    // the cylinder has its center of the base at x = 1.5, y = 2.5, z = 3.5 (in nm)
    // the cylinder has radius of 2.1 nm and height of 4.3 nm
    // the cylinder has its main axis oriented along the z-axis of the system
    let cylinder = Cylinder::new([1.5, 2.5, 3.5].into(), 2.1, 4.3, Dimension::Z);

    for atom in system.atoms_iter() {
        if cylinder.inside(atom.get_position().unwrap(), system.get_box_as_ref().unwrap()) {
            inside_cylinder.push(atom.clone());
        }
    }

    // atoms in the `inside_cylinder` vector are fully independent copies
    // of the original atoms in the `System` structure

    // print the number of extracted atoms
    println!("{:?}", inside_cylinder.len());

    Ok(())
}

函数式方法

use groan_rs::prelude::*;
use std::error::Error;

fn main() -> Result<(), Box<dyn Error + Send + Sync>> {
    // read a gro file
    let system = System::from_file("system.gro")?;

    // extract atoms that are located inside a specified cylinder
    let inside_cylinder: Vec<Atom> = system
        .atoms_iter()
        // the cylinder has its center of the base at x = 1.5, y = 2.5, z = 3.5 (in nm)
        // the cylinder has radius of 2.1 nm and height of 4.3 nm
        // the cylinder has its main axis oriented along the z-axis of the system
        .filter_geometry(Cylinder::new(
            [1.5, 2.5, 3.5].into(),
            2.1,
            4.3,
            Dimension::Z,
        ))
        .cloned()
        .collect();

    // atoms in the `inside_cylinder` vector are fully independent copies
    // of the original atoms in the `System` structure

    // print the number of extracted atoms
    println!("{}", inside_cylinder.len());

    Ok(())
}

groan选择语言

groan选择语言（GSL）是用于在groan_rs库中选择原子的查询语言。一般来说，GSL与VMD使用的选择语言类似。

基本查询

您可以根据以下条件选择原子：

• 使用resname XYZ根据它们的残基名称。例如，resname POPE将选择所有命名为POPE的系统中的原子。
• 使用resid XYZ或resnum XYZ根据它们的残基编号。例如，resid 17将选择所有编号为17的系统中的原子。
• 使用name XYZ或atomname XYZ根据它们的原子名称。例如，name P将选择所有名称为P的系统中的原子。
• 使用serial XYZ根据它们的真实原子编号。例如，serial 256将选择编号为256的原子（这保证是单个原子）。请注意，在这种情况下，原子是根据gromacs理解的“真实”原子编号选择的，而不是系统起源于的gro或pdb文件中指定的原子编号。
• 他们使用 atomid XYZ 或 atomnum XYZ 来指定原子序数。例如，atomid 124 将选择所有在系统来源的 gro 或 pdb 文件中原子序数为 124 的原子。这可能涉及多个原子。
• 他们使用 chain X 指定链标识符。例如 chain A 将选择属于链 'A' 的所有原子。链的信息通常存在于 pdb 文件中。注意，如果没有链信息，使用此关键字将不会选择任何原子。
• 他们使用 element name XYZ 或 elname XYZ 指定元素名称。例如，element name carbon 将选择所有碳原子。注意，原子必须被分配元素才能被选择，请参阅下文。
• 他们使用 element symbol X 或 elsymbol X 指定元素符号。例如，element symbol C 将选择所有碳原子。注意，原子必须被分配元素才能被选择，请参阅下文。
• 他们使用 label XYZ 指定标签，见下文。

use groan_rs::errors::WritePdbError;