Python写的调用Gromacs分子动力学MD的麦克斯韦速度分布测试类

代码：

class MaxwellBoltzmannTest(Test):
    @classmethod
    def parser(cls):
        parser = argparse.ArgumentParser(
            description='Tests the validity of the kinetic energy ensemble.',
            formatter_class=argparse.RawTextHelpFormatter,
            prog=cls.__name__
        )
        parser.add_argument('-t', '--tolerance', type=float, default=0.05,
                            help='The alpha value (confidence interval) used in\n'
                                 'the statistical test. Default: 0.05.')

        return parser

    @classmethod
    def prepare_parser(cls, input_dir, target_dir, system_name, nobackup, args):
        return cls.prepare(input_dir, target_dir, system_name, nobackup)

    @classmethod
    def analyze_parser(cls, gmx_parser, system_dir, system_name, base_data, verbosity, args):
        args = cls.parser().parse_args(args)
        return cls.analyze(gmx_parser, system_dir, system_name, base_data, verbosity,
                           alpha=args.tolerance)

    @classmethod
    def prepare(cls, input_dir, target_dir, system_name, nobackup):
        # no additional sims needed, base is enough
        # could check energy writing settings
        return []

    @classmethod
    def analyze(cls, gmx_parser, system_dir, system_name, base_data, verbosity, alpha=None):
        # Standard value
        if alpha is None:
            alpha = cls.parser().get_default('tolerance')   # alpha 是 tolerance

        # base data
        if base_data['reduced'] is None:
            current_dir = os.path.join(system_dir, 'base')
            base_data['reduced'] = gmx_parser.get_simulation_data(
                mdp=os.path.join(current_dir, 'mdout.mdp'),
                top=os.path.join(current_dir, 'system.top'),
                edr=os.path.join(current_dir, 'system.edr'),
                gro=os.path.join(current_dir, 'system.gro')
            )
        base_result = base_data['reduced']     # 关键数据

        p = kinetic_energy.mb_ensemble(base_result, verbosity=verbosity)  # 获取Maxwell速度概率
        # filename=os.path.join(system_dir, system_name + '_mb'))

        if p >= alpha:
            message = 'MaxwellBoltzmannTest PASSED (p = {:g}, alpha = {:f})'.format(p, alpha)
        else:
            message = 'MaxwellBoltzmannTest FAILED (p = {:g}, alpha = {:f})'.format(p, alpha)

        return {'test': p >= alpha,
                'result': p,
                'tolerance': alpha,
                'message': message}