哲●科●文Jerkwin 类别 标签 链接 关于 订阅

GROMACS编写轨迹分析工具的示例代码(gmx2018.2)

类别: 科   标签: gmx   阅读次数:   版权: (CC) BY-NC-SA ✕

GROMACS安装包括了一个使用轨迹分析框架编写轨迹分析工具的模板。

它可以从安装目录下的share/gromacs/template/和源代码分发中的share/template/中找到。

该文档的完整源代码也包含在本文档中：template.cpp本页的其余部分将介绍代码以解释不同的部分。

全局定义

首先包括一些通用的C++头文件：

#include <string>
#include <vector>

并继续包括分析库的头文件：

#include <gromacs/trajectoryanalysis.h>

此头文件包含其他头文件，这些头文件一起定义了编写轨迹分析工具所需的所有基本数据类型。为方便起见，我们还将gmx命名空间中的所有名称导入全局范围，以避免在任何地方重复使用该名称：

using namespace gmx;

工具模块类声明

然后我们定义一个实现分析工具的类：

class AnalysisTemplate : public TrajectoryAnalysisModule
{
	public:
		AnalysisTemplate();
		virtual void initOptions(IOptionsContainer          *options,
								 TrajectoryAnalysisSettings *settings);
		virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
								  const TopologyInformation        &top);
		virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
								  TrajectoryAnalysisModuleData *pdata);
		virtual void finishAnalysis(int nframes);
		virtual void writeOutput();
	private:
		class ModuleData;
		std::string                      fnDist_;
		double                           cutoff_;
		Selection                        refsel_;
		SelectionList                    sel_;
		AnalysisNeighborhood             nb_;
		AnalysisData                     data_;
		AnalysisDataAverageModulePointer avem_;
};

分析工具类来自于gmx::TrajectoryAnalysisModule，这个接口具有一些便捷的函数用于与其他代码早期接口进行对接.下面，我们将介绍模板中实现的不同方法（请注意，模板不会实现某些虚拟方法，因为它们不太常用），讨论在更复杂的情况下可能出现的一些问题。 有关可用虚拟方法和便捷功能的完整说明，请参阅gmx::TrajectoryAnalysisModule的文档。第一个成员变量块用于包含提供给不同选项的值。它们将根据分析工具的需要而变化。AnalysisNeighborhood对象提供在分析中使用的邻域搜索。最后一个变量块用于处理输出数据。有关如何使用它们的详细信息，请参照initAnalysis()。

对于模板，我们不需要任何自定义帧本地数据。 如果您认为需要一些更复杂的分析需求，请参阅gmx::TrajectoryAnalysisModuleData的文档以获取更多详细信息。 如果您不关心并行化，则无需考虑此部分。 您可以简单地在模块类中声明所有变量，在gmx::TrajectoryAnalysisModule::initAnalysis()中初始化它们，并在gmx::TrajectoryAnalysisModule::finishAnalysis()中进行任何后处理。

构造

分析模块的构造函数（和可能的析构函数）应该很简单：构造函数应该只是初始化默认值，析构函数应该释放模块管理的任何内存。对于模板，我们的类中没有需要显式释放的属性，因此我们只声明一个构造函数：

AnalysisTemplate::AnalysisTemplate()
: cutoff_(0.0)
{
registerAnalysisDataset(&data_, "avedist");
}

输入选项

模块的初始化分为几个方法，其中两个在模板中使用。gmx::TrajectoryAnalysisModule::initOptions()用于设置模块理解的选项，以及通过gmx::TrajectoryAnalysisSettings设置不同的选项（有关更多详细信息，请参阅该类的文档）：

void
AnalysisTemplate::initOptions(IOptionsContainer          *options,
							TrajectoryAnalysisSettings *settings)
{
	static const char *const desc[] = {
		"This is a template for writing your own analysis tools for",
		"GROMACS. The advantage of using GROMACS for this is that you",
		"have access to all information in the topology, and your",
		"program will be able to handle all types of coordinates and",
		"trajectory files supported by GROMACS. In addition,",
		"you get a lot of functionality for free from the trajectory",
		"analysis library, including support for flexible dynamic",
		"selections. Go ahead an try it![PAR]",
		"To get started with implementing your own analysis program,",
		"follow the instructions in the README file provided.",
		"This template implements a simple analysis programs that calculates",
		"average distances from a reference group to one or more",
		"analysis groups."
	};
	settings->setHelpText(desc);
	options->addOption(FileNameOption("o")
						.filetype(eftPlot).outputFile()
						.store(&fnDist_).defaultBasename("avedist")
						.description("Average distances from reference group"));
	options->addOption(SelectionOption("reference")
						.store(&refsel_).required()
						.description("Reference group to calculate distances from"));
	options->addOption(SelectionOption("select")
						.storeVector(&sel_).required().multiValue()
						.description("Groups to calculate distances to"));
	options->addOption(DoubleOption("cutoff").store(&cutoff_)
						.description("Cutoff for distance calculation (0 = no cutoff)"));
	settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
}

对于模板，我们首先为工具设置描述文本（用于帮助文本）。然后我们声明一个选项来指定输出文件名，后跟用于设置选择的选项，最后是一个设置截止值的选项。对于cutoff，默认值将是在构造函数中设置的值，但也可以在此处显式设置它。用户为选项提供的值将存储在成员变量中。最后，我们指出该工具始终需要拓扑信息。这样做仅用于演示目的：即使没有拓扑，模板中的代码也能正常工作。

有关如何定义不同类型选项的其他文档，请参阅gmx::IOptionsContainer，basicoptions.h和gmx::SelectionOption。 您只需要定义特定于分析的选项; 通用选项，例如，用于指定输入拓扑和添加轨迹框架。

要根据选项值调整设置或选择选项（例如，接受的选择数），您需要覆盖gmx::TrajectoryAnalysisModule::optionsFinished()。 为简单起见，这不是在模板中完成的。

分析初始化

实际分析在gmx::TrajectoryAnalysisModule::initAnalysis()中初始化：

void
AnalysisTemplate::initAnalysis(const TrajectoryAnalysisSettings &settings,
							   const TopologyInformation         & /*top*/)
{
	nb_.setCutoff(cutoff_);
	data_.setColumnCount(0, sel_.size());
	avem_.reset(new AnalysisDataAverageModule());
	data_.addModule(avem_);
	if (!fnDist_.empty())
	{
		AnalysisDataPlotModulePointer plotm(
				new AnalysisDataPlotModule(settings.plotSettings()));
		plotm->setFileName(fnDist_);
		plotm->setTitle("Average distance");
		plotm->setXAxisIsTime();
		plotm->setYLabel("Distance (nm)");
		data_.addModule(plotm);
	}
}

有关拓扑的信息作为参数传递。设置对象还可用于访问有关用户输入的信息。

此方法的主要任务之一是为它们设置适当的gmx::AnalysisData对象和模块（有关常规方法，请参阅gmx::TrajectoryAnalysisModule）。这些对象将用于处理工具的输出。 它们的主要目的是支持并行化，但即使您不关心并行性，它们仍然提供方便的构建块，例如，用于直方图和文件输出。

对于模板，我们首先设置邻域搜索的cutoff。

然后，我们创建并注册一个gmx::AnalysisData对象，该对象将为每个帧包含每个输入选择的一列。这将包含工具的主要输出：参考选择与特定选择之间的最小距离。然后，我们创建并设置一个模块，该模块将计算每个选择的平均距离（请参阅writeOutput()了解其使用方式）。最后，如果提供了输出文件，我们将创建并设置一个模块，该模块将绘制到文件的每帧距离。

如果分析模块在处理帧期间需要一些临时存储（即，它使用从gmx::TrajectoryAnalysisModuleData派生的自定义类），如果并行化被支持则应该在gmx::TrajectoryAnalysisModule::startFrames()（见下文）中分配。

如果需要根据第一帧中的数据进行初始化（最常见的是，基于框大小），则需要覆盖gmx::TrajectoryAnalysisModule::initAfterFirstFrame()，但这不在模板中使用。

分析框架

还有一个初始化方法需要重写以支持自动并行化：gmx::TrajectoryAnalysisModule::startFrames()。如果您不需要自定义帧本地数据（或根本不需要并行化），则可以跳过此方法并忽略gmx::TrajectoryAnalysisModule::analyzeFrame()的最后一个参数，以使事情更简单。在模板中，此方法不是必需的。

分析的主要部分是（在大多数分析代码中）在gmx::TrajectoryAnalysisModule::analyzeFrame()方法中完成，每个帧调用一次：

void
AnalysisTemplate::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
{

frnr参数给出了当前帧的从零开始的索引（主要用于gmx::AnalysisData），pbc包含用于距离计算的当前帧的PBC信息，例如，pbc_dx()，并且pdata指向在startFrames()中创建的数据结构。虽然通常不需要（时间字段除外），但可以通过fr访问原始帧数据。在大多数情况下，应该编写分析，使其通过选择获得所有位置数据，并且不假设它们具有恒定的大小。这就是支持选择引擎的完全灵活性所需的全部内容。

对于模板，我们首先从我们的自定义数据结构中获取数据以进行速记访问（如果使用自定义数据对象，则需要static_cast）：

AnalysisDataHandle     dh = pdata->dataHandle(data_);
const Selection        &refsel = pdata->parallelSelection(refsel_);

然后我们进行简单的计算并使用AnalysisDataHandle类来设置工具的每帧输出：

AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, refsel);
dh.startFrame(frnr, fr.time);
for (size_t g = 0; g < sel_.size(); ++g)
{
	const Selection &sel   = pdata->parallelSelection(sel_[g]);
	int              nr    = sel.posCount();
	real             frave = 0.0;
	for (int i = 0; i < nr; ++i)
	{
		SelectionPosition p = sel.position(i);
		frave += nbsearch.minimumDistance(p.x());
	}
	frave /= nr;
	dh.setPoint(g, frave);
}
dh.finishFrame();

处理完所有帧后，调用gmx::TrajectoryAnalysisModule::finishAnalysis()一次。这是对数据进行任何自定义后处理的地方。对于模板，我们什么都不做，因为所有必要的处理都在数据模块中完成：

void
AnalysisTemplate::finishAnalysis(int /*nframes*/)
{
}

如果在gmx::TrajectoryAnalysisModule::startFrames()中创建的数据结构用于跨帧聚合数据，则需要覆盖gmx::TrajectoryAnalysisModule::finishFrames()以组合数据结构中的数据（详情请参阅方法文档）。这对于模板来说不是必需的，因为ModuleData结构仅包含在分析单个帧期间使用的数据。

输出

最后，大多数程序需要在分析完成后输出一些值。在某些情况下，这可以通过正确链接数据模块来实现，但通常需要进行一些自定义处理。所有这些活动都应该在gmx::TrajectoryAnalysisModule::writeOutput()中完成。这使得更容易在例如脚本语言中重用分析模块，其中可能不期望输出到文件中。模板只打印出每个分析组的平均距离：

void
AnalysisTemplate::writeOutput()
{
	// We print out the average of the mean distances for each group.
	for (size_t g = 0; g < sel_.size(); ++g)
	{
		fprintf(stderr, "Average mean distance for '%s': %.3f nm\n",
				sel_[g].name(), avem_->average(0, g));
	}
}

在这里，我们使用avem_模块，我们在initAnalysis()中初始化它来聚合计算距离的平均值。

main()的定义

现在，唯一剩下的就是定义main()函数。 要实现命令行工具，它应该创建一个模块并使用gmx::TrajectoryAnalysisCommandLineRunner使用下面的样板代码运行它：

int
main(int argc, char *argv[])
{
	return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<AnalysisTemplate>(argc, argv);
}

GROMACS内置工具

使用模板实现的分析工具也可以轻松地包含在GROMACS库中。为此，请按照下列步骤操作：

1. 将您的工具源代码放入src/gromacs/trajectoryanalysis/modules/。

2. 删除using namespace gmx;, 并将所有代码包含在gmx::analysismodules命名空间中，并将工具类放入其中未命名的命名空间中。

3. 创建与您的工具对应的头文件，并使用gmx::analysismodules命名空间将以下类添加到其中（将Template替换为您的工具名称）：

class TemplateInfo
{
	public:
		static const char name[];
		static const char shortDescription[];
		static TrajectoryAnalysisModulePointer create();
};

4. 在源文件中，在未命名的命名空间之外添加这些项的定义（替换Template，AnalysisTemplate和具有正确值的字符串）：

const char TemplateInfo::name[]             = "template";
const char TemplateInfo::shortDescription[] =
	"Compute something";
TrajectoryAnalysisModulePointer TemplateInfo::create()
{
	return TrajectoryAnalysisModulePointer(new AnalysisTemplate);
}

5. 在src/gromacs/trajectoryanalysis/modules.cpp中注册您的模块。

6. 完成。现在可以使用您指定的名称作为gmx模板调用您的工具。

有关具体示例和文件的首选布局，请参阅src/gromacs/trajectoryanalysis/modules/中的现有工具。请使用现有文件中的Doxygen评论将自己记录为文件的作者。

◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆
前一篇:　GROMACS轨迹分析框架(gmx2018.2)
后一篇:　GROMACS团簇分析代码
访问人次(2015年7月 9日起): | 最后更新: 2024-04-16 06:38:20 UTC | 版权所有 © 2008 - 2024 Jerkwin 个