tests/physicalvalidation/gmx_physicalvalidation.py

   1 from __future__ import print_function, division, absolute_import
   2
   3
   4 import sys
   5 import os
   6 import shutil
   7 import json
   8 import argparse
   9 import re
  10 import math
  11 from collections import OrderedDict
  12
  13 from physical_validation import integrator, ensemble, kinetic_energy
  14 from physical_validation.util.gromacs_interface import GromacsInterface
  15 from physical_validation.data.gromacs_parser import GromacsParser
  16
  17
  18 def mkdir_bk(dirname, verbose=False, nobackup=False):
  19     if os.path.exists(dirname) and nobackup:
  20         shutil.rmtree(dirname)
  21     elif os.path.exists(dirname):
  22         if verbose:
  23             print('Directory ' + dirname + ' exists. Backing it up.')
  24         basename = os.path.basename(dirname)
  25         n = 1
  26         bk_dir = '#' + basename + '_' + str(n) + '#'
  27         while os.path.exists(dirname.replace(basename, bk_dir)):
  28             n += 1
  29             bk_dir = '#' + basename + '_' + str(n) + '#'
  30         os.rename(dirname, dirname.replace(basename, bk_dir))
  31     os.makedirs(dirname)
  32
  33
  34 def file_bk(filename, verbose=False):
  35     if os.path.exists(filename):
  36         if verbose:
  37             print('File ' + filename + ' exists. Backing it up.')
  38         basename = os.path.basename(filename)
  39         n = 1
  40         bk_file = '#' + basename + '_' + str(n) + '#'
  41         while os.path.exists(filename.replace(basename, bk_file)):
  42             n += 1
  43             bk_file = '#' + basename + '_' + str(n) + '#'
  44         os.rename(filename, filename.replace(basename, bk_file))
  45
  46
  47 def basic_run_cmds(directory, grompp_args=None, mdrun_args=None):
  48     grompp = '$GROMPPCMD -f system.mdp -p system.top -c system.gro -o system.tpr'
  49     if grompp_args:
  50         for arg in grompp_args:
  51             grompp += ' ' + arg
  52     mdrun = '$MDRUNCMD -s system.tpr -deffnm system'
  53     if mdrun_args:
  54         for arg in mdrun_args:
  55             mdrun += ' ' + arg
  56     return [
  57         'oldpath=$PWD',
  58         'cd ' + directory,
  59         grompp + ' && ' + mdrun,
  60         'cd $oldpath'
  61     ]
  62
  63
  64 class Test(object):
  65     @classmethod
  66     def parser(cls):
  67         raise NotImplementedError
  68
  69     @classmethod
  70     def prepare_parser(cls, input_dir, target_dir, system_name, nobackup, args):
  71         raise NotImplementedError
  72
  73     @classmethod
  74     def analyze_parser(cls, gmx_parser, system_dir, system_name, base_data, verbosity, args):
  75         raise NotImplementedError
  76
  77     @classmethod
  78     def prepare(cls, input_dir, target_dir, system_name, nobackup):
  79         raise NotImplementedError
  80
  81     @classmethod
  82     def analyze(cls, gmx_parser, system_dir, system_name, base_data, verbosity):
  83         raise NotImplementedError
  84
  85
  86 class IntegratorTest(Test):
  87     @classmethod
  88     def parser(cls):
  89         parser = argparse.ArgumentParser(
  90             description='Tests the integrator convergence.',
  91             formatter_class=argparse.RawTextHelpFormatter,
  92             prog=cls.__name__
  93         )
  94         parser.add_argument('-n', '--n_iterations', type=int, default=2,
  95                             help='The number of different time steps tested.')
  96         parser.add_argument('-t', '--tolerance', type=float, default=0.1,
  97                             help=('The relative tolerance accepted. Default: 0.1.\n'
  98                                   'Example: By default, the test passes if\n'
  99                                   '         3.6 <= dE(2*dt)/dE(dt) <= 4.4.')
 100                             )
 101
 102         return parser
 103
 104     @classmethod
 105     def prepare_parser(cls, input_dir, target_dir, system_name, nobackup, args):
 106         args = cls.parser().parse_args(args)
 107         return cls.prepare(input_dir, target_dir, system_name, nobackup,
 108                            n_iterations=args.n_iterations)
 109
 110     @classmethod
 111     def analyze_parser(cls, gmx_parser, system_dir, system_name, base_data, verbosity, args):
 112         args = cls.parser().parse_args(args)
 113         return cls.analyze(gmx_parser, system_dir, system_name, base_data, verbosity,
 114                            tolerance=args.tolerance, n_iterations=args.n_iterations)
 115
 116     @classmethod
 117     def prepare(cls, input_dir, target_dir, system_name, nobackup, n_iterations=None):
 118         # Standard value
 119         if n_iterations is None:
 120             n_iterations = cls.parser().get_default('n_iterations')
 121         # Read base options
 122         options = GromacsInterface.read_mdp(os.path.join(input_dir, 'system.mdp'))
 123
 124         # Check if options relevant to integrator tests are set
 125         # Set to standard if not - I am not happy about hardcoding this here...
 126         # For future: Find a way to obtain standards from GROMACS build
 127         if 'nsteps' not in options:
 128             raise ValueError('System ' + system_name + ' has no \'nsteps\' defined in mdp file. ' +
 129                              'Running integrator test does not make sense.')
 130         if 'dt' not in options:
 131             options['dt'] = str(0.001)
 132         if 'nstcalcenergy' not in options:
 133             options['nstcalcenergy'] = str(100)
 134         if 'nstenergy' not in options:
 135             options['nstenergy'] = str(1000)
 136         # if 'nstlist' not in options:
 137         #     options['nstlist'] = str(10)
 138
 139         # Prepare folders for iterations
 140         directories = []
 141         for n in range(1, n_iterations+1):
 142             current_dir = os.path.join(target_dir, 'integrator_' + str(n))
 143             mkdir_bk(current_dir, nobackup=nobackup)
 144             # update timesteps, length and intervals
 145             options['dt'] = str(float(options['dt'])*0.5)
 146             for key in ['nsteps', 'nstcalcenergy', 'nstenergy']:  # , 'nstlist']:
 147                 options[key] = str(int(options[key])*2)
 148             # write mdp
 149             GromacsInterface.write_mdp(options, os.path.join(current_dir, 'system.mdp'))
 150             # copy topology and starting structure
 151             for suffix in ['gro', 'top']:
 152                 shutil.copy2(os.path.join(input_dir, 'system.' + suffix),
 153                              current_dir)
 154             directories.append(current_dir)
 155
 156         return directories
 157
 158     @classmethod
 159     def analyze(cls, gmx_parser, system_dir, system_name, base_data, verbosity,
 160                 tolerance=None, n_iterations=None):
 161         # Standard value
 162         if tolerance is None:
 163             tolerance = cls.parser().get_default('tolerance')
 164         if n_iterations is None:
 165             n_iterations = cls.parser().get_default('n_iterations')
 166
 167         # list of results from simulations at different time steps
 168         results = []
 169
 170         # base data
 171         if base_data['reduced'] is None:
 172             current_dir = os.path.join(system_dir, 'base')
 173             base_data['reduced'] = gmx_parser.get_simulation_data(
 174                 mdp=os.path.join(current_dir, 'mdout.mdp'),
 175                 top=os.path.join(current_dir, 'system.top'),
 176                 edr=os.path.join(current_dir, 'system.edr'),
 177                 gro=os.path.join(current_dir, 'system.gro')
 178             )
 179         results.append(base_data['reduced'])
 180
 181         # append data at reduced dt
 182         for n in range(1, n_iterations+1):
 183             current_dir = os.path.join(system_dir, 'integrator_' + str(n))
 184             results.append(gmx_parser.get_simulation_data(
 185                 mdp=os.path.join(current_dir, 'mdout.mdp'),
 186                 top=os.path.join(current_dir, 'system.top'),
 187                 edr=os.path.join(current_dir, 'system.edr'),
 188                 gro=os.path.join(current_dir, 'system.gro')
 189             ))
 190
 191         # run test
 192         deviation = integrator.convergence(simulations=results,
 193                                            verbose=(verbosity > 0),
 194                                            convergence_test='max_deviation')
 195
 196         result = deviation <= tolerance
 197         if result:
 198             message = 'IntegratorTest PASSED (max deviation = {:.2g}, tolerance = {:.2g})'.format(
 199                 deviation, tolerance
 200             )
 201         else:
 202             message = 'IntegratorTest FAILED (max deviation = {:.2g}, tolerance = {:.2g})'.format(
 203                 deviation, tolerance
 204             )
 205
 206         return {'test': result,
 207                 'result': deviation,
 208                 'tolerance': tolerance,
 209                 'message': message}
 210
 211
 212 class EnsembleTest(Test):
 213     @classmethod
 214     def parser(cls):
 215         parser = argparse.ArgumentParser(
 216             description='Tests the validity of the potential energy and / or volume ensemble.',
 217             formatter_class=argparse.RawTextHelpFormatter,
 218             prog=cls.__name__
 219         )
 220         parser.add_argument('--dtemp', nargs='*', type=float, default=None,
 221                             help='Ensemble validations are made between two simulations at\n'
 222                                  'different state points.\n'
 223                                  'dtemp determines the temperature difference between base\n'
 224                                  'simulation and the additional point. If more than one\n'
 225                                  'value is given, several tests will be performed.\n'
 226                                  'By also giving dpress, both temperature and pressure can\n'
 227                                  'be displaced simultaneously.')
 228         parser.add_argument('--dpress', nargs='*', type=float, default=None,
 229                             help='Ensemble validations are made between two simulations at\n'
 230                                  'different state points.\n'
 231                                  'dpress determines the pressure difference between base\n'
 232                                  'simulation and the additional point. If more than one\n'
 233                                  'value is given, several tests will be performed.\n'
 234                                  'By also giving dtemp, both temperature and pressure can\n'
 235                                  'be displaced simultaneously.')
 236         parser.add_argument('-t', '--tolerance', type=float, default=3,
 237                             help=('The number of standard deviations a result can be off\n'
 238                                   'to be still accepted. Default: 3.'))
 239
 240         return parser
 241
 242     @classmethod
 243     def prepare_parser(cls, input_dir, target_dir, system_name, nobackup, args):
 244         args = cls.parser().parse_args(args)
 245         return cls.prepare(input_dir, target_dir, system_name, nobackup,
 246                            dtemp=args.dtemp, dpress=args.dpress)
 247
 248     @classmethod
 249     def analyze_parser(cls, gmx_parser, system_dir, system_name, base_data, verbosity, args):
 250         args = cls.parser().parse_args(args)
 251         return cls.analyze(gmx_parser, system_dir, system_name, base_data, verbosity,
 252                            tolerance=args.tolerance, dtemp=args.dtemp, dpress=args.dpress)
 253
 254     @classmethod
 255     def prepare(cls, input_dir, target_dir, system_name, nobackup, dtemp=None, dpress=None):
 256         # No standard values (system-dependent!)
 257         if not dtemp and not dpress:
 258             raise ValueError('Ensemble test for system ' + system_name +
 259                              ' has no defined temperature or pressure difference.')
 260         # Pad arrays - assume 0 difference if other difference is set
 261         if not dtemp:
 262             dtemp = [0] * len(dpress)
 263         if not dpress:
 264             dpress = [0] * len(dtemp)
 265         if len(dtemp) < len(dpress):
 266             dtemp.extend([0] * (len(dpress) - len(dtemp)))
 267         if len(dpress) < len(dtemp):
 268             dpress.extend([0] * (len(dtemp) - len(dpress)))
 269         # Check if we do any pressure differences
 270         no_press = True
 271         for dp in dpress:
 272             if dp*dp >= 1e-12:
 273                 no_press = False
 274                 break
 275
 276         # Read base options
 277         options = GromacsInterface.read_mdp(os.path.join(input_dir, 'system.mdp'))
 278
 279         # Check for target temperature
 280         if 'ref-t' in options:
 281             ref_t = float(options['ref-t'])
 282             ref_t_key = 'ref-t'
 283         elif 'ref_t' in options:
 284             ref_t = float(options['ref_t'])
 285             ref_t_key = 'ref_t'
 286         else:
 287             raise ValueError('Ensemble test for system ' + system_name + ' selected, ' +
 288                              'but system has no defined temperature.')
 289
 290         # Check for target pressure, if relevant
 291         if no_press:
 292             ref_p = None
 293             ref_p_key = None
 294         elif 'ref-p' in options:
 295             ref_p = float(options['ref-p'])
 296             ref_p_key = 'ref-p'
 297         elif 'ref_p' in options:
 298             ref_p = float(options['ref_p'])
 299             ref_p_key = 'ref_p'
 300         else:
 301             raise ValueError('Ensemble test with pressure difference for system ' + system_name +
 302                              ' selected, but system has no defined pressure.')
 303
 304         # Loop over tests
 305         directories = []
 306         for n, (dt, dp) in enumerate(zip(dtemp, dpress)):
 307             current_dir = os.path.join(target_dir, 'ensemble_' + str(n+1))
 308             mkdir_bk(current_dir, nobackup=nobackup)
 309             # change temperature & pressure
 310             options[ref_t_key] = str(ref_t + dt)
 311             if not no_press:
 312                 options[ref_p_key] = str(ref_p + dp)
 313             # write mdp
 314             GromacsInterface.write_mdp(options, os.path.join(current_dir, 'system.mdp'))
 315             # copy topology and starting structure
 316             for suffix in ['gro', 'top']:
 317                 shutil.copy2(os.path.join(input_dir, 'system.' + suffix),
 318                              current_dir)
 319             directories.append(current_dir)
 320
 321         return directories
 322
 323     @classmethod
 324     def analyze(cls, gmx_parser, system_dir, system_name, base_data, verbosity,
 325                 tolerance=None, dtemp=None, dpress=None):
 326         # No standard values (system-dependent!)
 327         if not dtemp and not dpress:
 328             raise ValueError('Ensemble test for system ' + system_name +
 329                              ' has no defined temperature or pressure difference.')
 330         # Pad arrays - assume 0 difference if other difference is set
 331         if not dtemp:
 332             dtemp = [0] * len(dpress)
 333         if not dpress:
 334             dpress = [0] * len(dtemp)
 335         if len(dtemp) < len(dpress):
 336             dtemp.extend([0] * (len(dpress) - len(dtemp)))
 337         if len(dpress) < len(dtemp):
 338             dpress.extend([0] * (len(dtemp) - len(dpress)))
 339
 340         nsystems = len(dtemp)
 341
 342         if tolerance is None:
 343             tolerance = cls.parser().get_default('tolerance')
 344
 345         # base data
 346         if base_data['reduced'] is None:
 347             current_dir = os.path.join(system_dir, 'base')
 348             base_data['reduced'] = gmx_parser.get_simulation_data(
 349                 mdp=os.path.join(current_dir, 'mdout.mdp'),
 350                 top=os.path.join(current_dir, 'system.top'),
 351                 edr=os.path.join(current_dir, 'system.edr'),
 352                 gro=os.path.join(current_dir, 'system.gro')
 353             )
 354         base_result = base_data['reduced']
 355
 356         # list of results from simulations at different state points
 357         results = []
 358         for n in range(nsystems):
 359             current_dir = os.path.join(system_dir, 'ensemble_' + str(n+1))
 360             results.append(gmx_parser.get_simulation_data(
 361                 mdp=os.path.join(current_dir, 'mdout.mdp'),
 362                 top=os.path.join(current_dir, 'system.top'),
 363                 edr=os.path.join(current_dir, 'system.edr'),
 364                 gro=os.path.join(current_dir, 'system.gro')
 365             ))
 366
 367         # run tests
 368         passed = True
 369         message = ''
 370         max_quantiles = -1
 371         for result, dt, dp in zip(results, dtemp, dpress):
 372             quantiles = ensemble.check(base_result, result, verbosity=verbosity)
 373             # filename=os.path.join(system_dir, system_name + '_ens'))
 374             if any(q > tolerance or math.isnan(q) for q in quantiles):
 375                 passed = False
 376                 if len(quantiles) == 1:
 377                     message += '\n    --dtemp={:.1f} --dpress={:.1f} : FAILED ({:.1f} quantiles off)'.format(
 378                         dt, dp, quantiles[0]
 379                     )
 380                 else:
 381                     message += '\n    --dtemp={:.1f} --dpress={:.1f} : FAILED ([{:.1f}, {:.1f}] quantiles off)'.format(
 382                         dt, dp, quantiles[0], quantiles[1]
 383                     )
 384             else:
 385                 if len(quantiles) == 1:
 386                     message += '\n    --dtemp={:.1f} --dpress={:.1f} : PASSED ({:.1f} quantiles off)'.format(
 387                         dt, dp, quantiles[0]
 388                     )
 389                 else:
 390                     message += '\n    --dtemp={:.1f} --dpress={:.1f} : PASSED ([{:.1f}, {:.1f}] quantiles off)'.format(
 391                         dt, dp, quantiles[0], quantiles[1]
 392                     )
 393             max_quantiles = max(max_quantiles, max(quantiles))
 394
 395         if passed:
 396             message = ('EnsembleTest PASSED (tolerance: {:.1f} quantiles)'.format(tolerance) +
 397                        message)
 398         else:
 399             message = ('EnsembleTest FAILED (tolerance: {:.1f} quantiles)'.format(tolerance) +
 400                        message)
 401
 402         return {'test': passed,
 403                 'result': max_quantiles,
 404                 'tolerance': tolerance,
 405                 'message': message}
 406
 407
 408 class MaxwellBoltzmannTest(Test):
 409     @classmethod
 410     def parser(cls):
 411         parser = argparse.ArgumentParser(
 412             description='Tests the validity of the kinetic energy ensemble.',
 413             formatter_class=argparse.RawTextHelpFormatter,
 414             prog=cls.__name__
 415         )
 416         parser.add_argument('-t', '--tolerance', type=float, default=0.05,
 417                             help='The alpha value (confidence interval) used in\n'
 418                                  'the statistical test. Default: 0.05.')
 419
 420         return parser
 421
 422     @classmethod
 423     def prepare_parser(cls, input_dir, target_dir, system_name, nobackup, args):
 424         return cls.prepare(input_dir, target_dir, system_name, nobackup)
 425
 426     @classmethod
 427     def analyze_parser(cls, gmx_parser, system_dir, system_name, base_data, verbosity, args):
 428         args = cls.parser().parse_args(args)
 429         return cls.analyze(gmx_parser, system_dir, system_name, base_data, verbosity,
 430                            alpha=args.tolerance)
 431
 432     @classmethod
 433     def prepare(cls, input_dir, target_dir, system_name, nobackup):
 434         # no additional sims needed, base is enough
 435         # could check energy writing settings
 436         return []
 437
 438     @classmethod
 439     def analyze(cls, gmx_parser, system_dir, system_name, base_data, verbosity, alpha=None):
 440         # Standard value
 441         if alpha is None:
 442             alpha = cls.parser().get_default('tolerance')
 443
 444         # base data
 445         if base_data['reduced'] is None:
 446             current_dir = os.path.join(system_dir, 'base')
 447             base_data['reduced'] = gmx_parser.get_simulation_data(
 448                 mdp=os.path.join(current_dir, 'mdout.mdp'),
 449                 top=os.path.join(current_dir, 'system.top'),
 450                 edr=os.path.join(current_dir, 'system.edr'),
 451                 gro=os.path.join(current_dir, 'system.gro')
 452             )
 453         base_result = base_data['reduced']
 454
 455         p = kinetic_energy.mb_ensemble(base_result, verbosity=verbosity)
 456         # filename=os.path.join(system_dir, system_name + '_mb'))
 457
 458         if p >= alpha:
 459             message = 'MaxwellBoltzmannTest PASSED (p = {:g}, alpha = {:f})'.format(p, alpha)
 460         else:
 461             message = 'MaxwellBoltzmannTest FAILED (p = {:g}, alpha = {:f})'.format(p, alpha)
 462
 463         return {'test': p >= alpha,
 464                 'result': p,
 465                 'tolerance': alpha,
 466                 'message': message}
 467
 468
 469 class EquipartitionTest(Test):
 470     @classmethod
 471     def parser(cls):
 472         parser = argparse.ArgumentParser(
 473             description='Tests the equipartition of the kinetic energy.',
 474             formatter_class=argparse.RawTextHelpFormatter,
 475             prog=cls.__name__
 476         )
 477         parser.add_argument('--distribution', default=False, action='store_true',
 478                             help=('If set, the groups of degrees of freedom will\n'
 479                                   'be separately tested as Maxwell-Boltzmann\n'
 480                                   'distributions. Otherwise, only the difference in\n'
 481                                   'average kinetic energy is checked.'))
 482         parser.add_argument('-t', '--tolerance', type=float, default=0.1,
 483                             help=('If --distribution is set:\n'
 484                                   '    The alpha value (confidence interval) used in the\n'
 485                                   '    statistical test. Default: 0.1.\n'
 486                                   'Else:\n'
 487                                   '    The maximum deviation allowed between groups.\n'
 488                                   '    Default: 0.1 (10%%)'))
 489         parser.add_argument('--random_groups', type=int, default=0,
 490                             help='Number of random division tests attempted.\n'
 491                                  'Default: 0 (random division tests off).')
 492         parser.add_argument('--random_divisions', type=int, default=2,
 493                             help='Number of groups the system is randomly divided in.\n'
 494                                  'Default: 2.')
 495         parser.add_argument('--molec_group', type=int, nargs='*', default=None, action='append',
 496                             help=('List of molecule indeces defining a group. Useful to\n'
 497                                   'pre-define groups of molecules\n'
 498                                   '(e.g. solute / solvent, liquid mixture species, ...).\n'
 499                                   'If not set, no pre-defined molecule groups will be\ntested.\n'
 500                                   'Note: If the last --molec-group flag is given empty,\n'
 501                                   'the remaining molecules are collected in this group.\n'
 502                                   'This allows, for example, to only specify the solute,\n'
 503                                   'and indicate the solvent by giving an empty flag.'))
 504
 505         return parser
 506
 507     @classmethod
 508     def prepare(cls, input_dir, target_dir, system_name, nobackup):
 509         # no additional sims needed, base is enough
 510         # could check position, velocity & energy writing settings
 511         return []
 512
 513     @classmethod
 514     def prepare_parser(cls, input_dir, target_dir, system_name, nobackup, args):
 515         return cls.prepare(input_dir, target_dir, system_name, nobackup)
 516
 517     @classmethod
 518     def analyze_parser(cls, gmx_parser, system_dir, system_name, base_data, verbosity, args):
 519         args = cls.parser().parse_args(args)
 520         return cls.analyze(gmx_parser, system_dir, system_name, base_data, verbosity,
 521                            tolerance=args.tolerance, distribution=args.distribution,
 522                            random_groups=args.random_groups, random_divisions=args.random_divisions,
 523                            molec_group=args.molec_group)
 524
 525     @classmethod
 526     def analyze(cls, gmx_parser, system_dir, system_name, base_data, verbosity,
 527                 tolerance=None, distribution=None,
 528                 random_groups=None, random_divisions=None,
 529                 molec_group=None):
 530
 531         # Standard values
 532         if tolerance is None:
 533             tolerance = cls.parser().get_default('tolerance')
 534         if distribution is None:
 535             distribution = cls.parser().get_default('distribution')
 536         if random_groups is None:
 537             random_groups = cls.parser().get_default('random_groups')
 538         if random_divisions is None:
 539             random_divisions = cls.parser().get_default('random_divisions')
 540         if molec_group is None:
 541             molec_group = cls.parser().get_default('molec_group')
 542
 543         # base data
 544         if base_data['full'] is None:
 545             current_dir = os.path.join(system_dir, 'base')
 546             base_data['full'] = gmx_parser.get_simulation_data(
 547                 mdp=os.path.join(current_dir, 'mdout.mdp'),
 548                 top=os.path.join(current_dir, 'system.top'),
 549                 edr=os.path.join(current_dir, 'system.edr'),
 550                 gro=os.path.join(current_dir, 'system.trr')
 551             )
 552         base_result = base_data['full']
 553
 554         result = kinetic_energy.equipartition(base_result,
 555                                               dtemp=tolerance, alpha=tolerance,
 556                                               distribution=distribution, molec_groups=molec_group,
 557                                               random_divisions=random_divisions, random_groups=random_groups,
 558                                               verbosity=0)
 559
 560         if distribution:
 561             res = min(result)
 562             test = res >= tolerance
 563             if test:
 564                 message = 'EquipartitionTest PASSED (min p = {:g}, alpha = {:f})'.format(res, tolerance)
 565             else:
 566                 message = 'EquipartitionTest FAILED (min p = {:g}, alpha = {:f})'.format(res, tolerance)
 567         else:
 568             dev_min = min(result)
 569             dev_max = max(result)
 570             if (1-dev_min) > (dev_max-1):
 571                 res = dev_min
 572             else:
 573                 res = dev_max
 574             test = res <= tolerance
 575             if test:
 576                 message = 'EquipartitionTest PASSED (max dev = {:g}, tolerance = {:f})'.format(res, tolerance)
 577             else:
 578                 message = 'EquipartitionTest FAILED (max dev = {:g}, tolerance = {:f})'.format(res, tolerance)
 579
 580         return {'test': test,
 581                 'result': res,
 582                 'tolerance': tolerance,
 583                 'message': message}
 584
 585
 586 class KinConstraintsTest(Test):
 587     @classmethod
 588     def parser(cls):
 589         parser = argparse.ArgumentParser(
 590             description='Tests whether there is kinetic energy in constrained degrees of\nfreedom.',
 591             formatter_class=argparse.RawTextHelpFormatter,
 592             prog=cls.__name__
 593         )
 594         parser.add_argument('-t', '--tolerance', type=float, default=0.1,
 595                             help='Tolerance - TODO. Default: 0.1.')
 596
 597         return parser
 598
 599     @classmethod
 600     def prepare(cls, input_dir, target_dir, system_name, nobackup):
 601         # no additional sims needed, base is enough
 602         # could check if there are any constraints in the system
 603         return []
 604
 605     @classmethod
 606     def prepare_parser(cls, input_dir, target_dir, system_name, nobackup, args):
 607         return cls.prepare(input_dir, target_dir, system_name, nobackup)
 608
 609     @classmethod
 610     def analyze_parser(cls, gmx_parser, system_dir, system_name, base_data, verbosity, args):
 611         raise NotImplementedError
 612
 613     @classmethod
 614     def analyze(cls, gmx_parser, system_dir, system_name, base_data, verbosity):
 615         raise NotImplementedError
 616
 617
 618 all_tests = OrderedDict([
 619     ('integrator', IntegratorTest),
 620     ('ensemble', EnsembleTest),
 621     ('kin_mb', MaxwellBoltzmannTest),
 622     ('kin_equipartition', EquipartitionTest),
 623     ('kin_constraints', KinConstraintsTest)
 624 ])
 625
 626
 627 def parse_systems(systems_json, systems_user, source_path,
 628                   analyze_only):
 629     # Parse json
 630     # As the order of the systems and the tests
 631     # might be meaningful, we need ordered dicts!
 632     system_list = json.load(systems_json)
 633     system_dict = OrderedDict()
 634     for system in system_list:
 635         system_name = system['name']
 636         system_dir = system['dir']
 637         # do the input files exist? (only relevant if we're not only analyzing)
 638         if not analyze_only:
 639             input_dir = os.path.join(source_path, system_dir, 'input')
 640             if not (os.path.isdir(input_dir) and
 641                     os.path.exists(os.path.join(input_dir, 'system.mdp')) and
 642                     os.path.exists(os.path.join(input_dir, 'system.top')) and
 643                     os.path.exists(os.path.join(input_dir, 'system.gro'))):
 644                 raise ValueError('System ' + system_name + ' in ' +
 645                                  systems_json.name + ': Input files not found')
 646         # no need to run systems that we don't test
 647         if 'tests' not in system:
 648             raise ValueError('System ' + system_name + ' in ' +
 649                              systems_json.name + ' has no tests defined')
 650         test_list = system['tests']
 651         system['tests'] = OrderedDict()
 652         for test in test_list:
 653             test_name = test['test']
 654             if test_name not in all_tests:
 655                 raise ValueError('Test ' + test_name + ' in ' +
 656                                  systems_json.name + ' is not a valid test')
 657             if test_name not in system['tests']:
 658                 if 'args' in test:
 659                     test['args'] = [test['args'].split()]
 660                 else:
 661                     test['args'] = [[]]
 662                 system['tests'][test_name] = test
 663             else:
 664                 if 'args' in test:
 665                     system['tests'][test_name]['args'].append(test['args'].split())
 666                 else:
 667                     system['tests'][test_name]['args'].append([])
 668
 669         system_dict[system_name] = system
 670         # add standard arguments
 671         if 'grompp_args' not in system:
 672             system['grompp_args'] = []
 673         else:
 674             system['grompp_args'] = system['grompp_args'].split()
 675         if 'mdrun_args' not in system:
 676             system['mdrun_args'] = []
 677         else:
 678             system['mdrun_args'] = system['mdrun_args'].split()
 679
 680     # if user has not chosen otherwise, return full dict of systems
 681     if not systems_user:
 682         return system_dict
 683
 684     # make sure user_dicts matches at least something
 685     for user_system in systems_user:
 686         for system in system_dict:
 687             if re.match(user_system + '$', system):
 688                 break
 689         else:
 690             raise ValueError('System ' + user_system +
 691                              ' used in command line argument is not defined in ' +
 692                              systems_json.name)
 693
 694     for system in list(system_dict.keys()):
 695         # delete systems not selected by user
 696         for user_system in systems_user:
 697             if re.match(user_system + '$', system):
 698                 break
 699         else:
 700             system_dict.pop(system)
 701
 702     # return reduced dict of systems
 703     return system_dict
 704
 705
 706 def main(args):
 707     parser = argparse.ArgumentParser(
 708         description='Physical validation suite for GROMACS.',
 709         prog='gmx_physicalvalidation.py',
 710         formatter_class=argparse.RawTextHelpFormatter,
 711         epilog='Use --tests for details about the available tests and their arguments.'
 712     )
 713     parser.add_argument('json', type=open,
 714                         metavar='systems.json',
 715                         help='Json file containing systems and tests to be ran.')
 716     parser.add_argument('--tests', default=False, action='store_true',
 717                         help='Print details about the available tests and their arguments and exit.')
 718     parser.add_argument('-v', '--verbosity', type=int,
 719                         metavar='v', default=0,
 720                         help='Verbosity level. Default: 0 (quiet).')
 721     parser.add_argument('--gmx', type=str, metavar='exe', default=None,
 722                         help=('GROMACS executable. Default: Trying to use \'gmx\'.\n' +
 723                               'Note: If -p is used, the GROMACS executable is not needed.'))
 724     parser.add_argument('--bindir', type=str, metavar='dir', default=None,
 725                         help=('GROMACS binary directory.\n' +
 726                               'If set, trying to use \'bindir/gmx\' instead of plain \'gmx\'\n' +
 727                               'Note: If --gmx is set, --bindir and --suffix are ignored.'))
 728     parser.add_argument('--suffix', type=str, metavar='_s', default=None,
 729                         help=('Suffix of the GROMACS executable.\n' +
 730                               'If set, trying to use \'gmx_s\' instead of plain \'gmx\'\n' +
 731                               'Note: If --gmx is set, --bindir and --suffix are ignored.'))
 732     group = parser.add_mutually_exclusive_group()
 733     group.add_argument('-p', '--prepare', action='store_true',
 734                        default=False,
 735                        help=('Only prepare simulations and output a \'run.sh\' file\n' +
 736                              'containing the necessary commands to run the systems.\n' +
 737                              'This allows to separate running simulations from analyzing them,\n' +
 738                              'useful e.g. to analyze the results on a different machine.\n' +
 739                              'Default: If none of \'-p\', \'-r\' or \'-a\' is given,\n' +
 740                              '         the systems are prepared, ran and analyzed in one call.'))
 741     group.add_argument('-r', '--run', action='store_true',
 742                        default=False,
 743                        help=('Only prepare and run simulations.\n' +
 744                              'Default: If none of \'-p\', \'-r\' or \'-a\' is given,\n' +
 745                              '         the systems are prepared, ran and analyzed in one call.'))
 746     group.add_argument('-a', '--analyze', action='store_true',
 747                        default=False,
 748                        help=('Only analyze previously ran simulations.\n' +
 749                              'This requires that the systems have been prepared using this program.\n' +
 750                              'Default: If none of \'-p\', \'-r\' or \'-a\' is given,\n' +
 751                              '         the systems are prepared, ran and analyzed in one call.'))
 752     parser.add_argument('-s', '--system', action='append', dest='systems',
 753                         metavar='system',
 754                         help=('Specify which system to run.\n' +
 755                               '\'system\' needs to match a system defined in the json file.\n' +
 756                               'Several systems can be specified by chaining several \'-s\' arguments.\n' +
 757                               'Defaults: If no \'-s\' argument is given, all systems and tests\n' +
 758                               '          defined in the json file are ran.\n' +
 759                               'Note: \'system\' can be a regular expression matching more than one system.'))
 760     parser.add_argument('--mpicmd', type=str, metavar='cmd', default=None,
 761                         help='MPI command used to invoke run command')
 762     parser.add_argument('--wd', '--working_dir', type=str,
 763                         metavar='dir', default=None,
 764                         help='Working directory (default: current directory)')
 765     parser.add_argument('--nobackup', default=False, action='store_true',
 766                         help='Do not create backups of files or folders.')
 767
 768     if '--tests' in args:
 769         message = ('Physical validation suite for GROMACS\n'
 770                    'Available tests and their arguments\n'
 771                    '=====================================\n\n'
 772                    'The following tests can be specified in the json file:\n\n')
 773         for test in all_tests:
 774             message += '  * ' + test + '\n'
 775
 776         message += '\nAll tests accept additional arguments:\n'
 777
 778         for test, test_cls in all_tests.items():
 779             message += '\n'
 780             test_help = test_cls.parser().format_help()
 781             test_help = test_help.replace(test_cls.__name__, test).replace('usage: ', '')
 782             test_help = test_help.replace('  -h, --help            show this help message and exit\n', '')
 783             test_help = test_help.replace(' [-h]', '')
 784             test_help = test_help.replace(' ' * (len(test_cls.__name__) + 7), ' ' * len(test))
 785
 786             first_line = test_help.split('\n')[0]
 787             separator = '-' * 70
 788             message += test_help.replace(first_line, separator + '\n' + first_line)
 789
 790         sys.stderr.write(message)
 791         return
 792
 793     args = parser.parse_args(args)
 794
 795     # the input files are expected to be located where this script is
 796     source_path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'systems')
 797     # target directory can be current or user chosen
 798     if args.wd is None:
 799         target_path = os.getcwd()
 800     else:
 801         if not os.path.exists(args.wd):
 802             os.makedirs(args.wd)
 803         target_path = args.wd
 804
 805     # parse simulation stage to perform
 806     do_all = not (args.prepare or args.run or args.analyze)
 807     do_prepare = do_all or args.prepare or args.run
 808     write_script = args.prepare
 809     do_run = do_all or args.run
 810     do_analysis = do_all or args.analyze
 811
 812     # get ordered dict of systems from combination of json file and user choices
 813     systems = parse_systems(args.json, args.systems, source_path, args.analyze)
 814
 815     # prepare GROMACS interface
 816     if args.gmx:
 817         gmx = args.gmx
 818     else:
 819         gmx = 'gmx'
 820         if args.suffix:
 821             gmx += args.suffix
 822         if args.bindir:
 823             gmx = os.path.join(args.bindir, gmx)
 824     gmx_interface = None
 825     gmx_parser = None
 826     if do_run or do_analysis:
 827         gmx_interface = GromacsInterface(exe=gmx)
 828         gmx_parser = GromacsParser(exe=gmx)
 829
 830     if do_prepare:
 831         nsystems = len(systems)
 832         n = 0
 833         runs = []  # this will contain all information needed to run the system
 834         for system_name, system in systems.items():
 835             n += 1
 836             print('\rPreparing run files for systems... [{:d}/{:d}] '.format(n, nsystems), end='')
 837             sys.stdout.flush()  # py2 compatibility
 838             system_dir = system['dir']
 839             system_dirs = []  # list of directories with subsystems
 840             # prepare the base system
 841             input_dir = os.path.join(source_path, system_dir, 'input')
 842             target_dir = os.path.join(target_path, system_dir)
 843             mkdir_bk(target_dir, nobackup=args.nobackup)
 844             basedir = os.path.join(target_dir, 'base')
 845             mkdir_bk(basedir, nobackup=args.nobackup)
 846             for suffix in ['mdp', 'gro', 'top']:
 847                 shutil.copy2(os.path.join(input_dir, 'system.' + suffix),
 848                              basedir)
 849             system_dirs.append(basedir)
 850             # call prepare method of chosen tests
 851             for test_name, test in system['tests'].items():
 852                 for test_args in test['args']:
 853                     system_dirs.extend(
 854                         all_tests[test_name].prepare_parser(input_dir, target_dir, system_name,
 855                                                             args.nobackup, test_args)
 856                     )
 857
 858             # save run information
 859             for d in system_dirs:
 860                 runs.append({
 861                     'dir': d,
 862                     'grompp_args': system['grompp_args'],
 863                     'mdrun_args': system['mdrun_args']
 864                 })
 865         # end of loop over systems
 866         print('-- done.')
 867
 868         if write_script:
 869             print('Writing run script... ', end='')
 870             sys.stdout.flush()  # py2 compatibility
 871             script_file = os.path.join(target_path, 'run_simulations.sh')
 872             if not args.nobackup:
 873                 file_bk(script_file)
 874             with open(script_file, 'w') as f:
 875                 f.write('# This file was created by the physical validation suite for GROMACS.\n')
 876                 f.write('\n# Define run variables\n')
 877                 f.write('WORKDIR=' + os.path.abspath(target_path) + '\n')
 878                 f.write('GROMPPCMD="' + os.path.abspath(gmx) + ' grompp"\n')
 879                 f.write('MDRUNCMD="' + os.path.abspath(gmx) + ' mdrun"\n')
 880                 f.write('\n# Run systems\n')
 881                 f.write('startpath=$PWD\n')
 882                 f.write('cd $WORKDIR\n')
 883                 for run in runs:
 884                     for cmd in basic_run_cmds(directory=os.path.relpath(os.path.abspath(run['dir']),
 885                                                                         os.path.abspath(target_path)),
 886                                               grompp_args=run['grompp_args'],
 887                                               mdrun_args=run['mdrun_args']):
 888                         f.write(cmd + '\n')
 889                     f.write('\n')
 890                 f.write('cd $startpath\n')
 891             print('-- done.')
 892             print('Run script written to ' + script_file)
 893             print('Adapt script as necessary and run simulations. Make sure to preserve the folder structure!')
 894             print('Once all simulations have ran, analyze the results using `make check-phys-analyze` or '
 895                   'using the `-a` flag of `gmx_physicalvalidation.py`.')
 896         # end if write_script
 897
 898         if do_run:
 899             nruns = len(runs)
 900             # send messages from GROMACS to log
 901             gmx_log = open(os.path.join(target_path, 'physicalvalidation_gmx.log'), 'w')
 902             for n, run in enumerate(runs):
 903                 print('\rRunning (sub)systems... [{:d}/{:d}] '.format(n+1, nruns), end='')
 904                 sys.stdout.flush()  # py2 compatibility
 905                 gmx_interface.grompp(mdp='system.mdp',
 906                                      top='system.top',
 907                                      gro='system.gro',
 908                                      tpr='system.tpr',
 909                                      cwd=run['dir'],
 910                                      args=run['grompp_args'],
 911                                      stdout=gmx_log,
 912                                      stderr=gmx_log)
 913                 gmx_interface.mdrun(tpr='system.tpr',
 914                                     deffnm='system',
 915                                     cwd=run['dir'],
 916                                     args=run['mdrun_args'],
 917                                     stdout=gmx_log,
 918                                     stderr=gmx_log,
 919                                     mpicmd=args.mpicmd)
 920             gmx_log.close()
 921             print('-- done.')
 922         # end if do_run
 923     # end if do_prepare
 924
 925     if do_analysis:
 926         title = 'GROMACS PHYSICAL VALIDATION RESULTS'
 927         width = 70
 928         indent = int((width - len(title))/2)
 929         print()
 930         print(' ' * indent + '=' * len(title))
 931         print(' ' * indent + title)
 932         print(' ' * indent + '=' * len(title))
 933         print()
 934         passed = True
 935         for system_name, system in systems.items():
 936             system_dir = system['dir']
 937             # save system data if re-used for different test
 938             # massively reduces run time of multiple tests
 939             system_data = {
 940                 'reduced': None,
 941                 'full': None
 942             }
 943             # system directory
 944             target_dir = os.path.join(target_path, system_dir)
 945
 946             print('Analyzing system ' + system_name)
 947
 948             # call analyze method of chosen tests
 949             for test_name, test in system['tests'].items():
 950                 for test_args in test['args']:
 951                     try:
 952                         result = all_tests[test_name].analyze_parser(gmx_parser, target_dir,
 953                                                                      system_name, system_data,
 954                                                                      args.verbosity, test_args)
 955                     except Exception as err:
 956                         print('    ' + all_tests[test_name].__name__ + ' FAILED (Exception in evaluation)')
 957                         print('    '*2 + type(err).__name__ + ': ' + str(err))
 958                         passed = False
 959                     else:
 960                         for line in result['message'].split('\n'):
 961                             print('    ' + line)
 962
 963                         passed = passed and result['test']
 964             # end loop over tests
 965             print()
 966         # end loop over systems
 967         return int(not passed)
 968     # end if do_analysis
 969
 970     # assuming everything is ok if we ended up here
 971     return 0
 972
 973
 974 if __name__ == "__main__":
 975     return_value = main(sys.argv[1:])
 976     sys.exit(return_value)