tests/physicalvalidation/gmx_physicalvalidation.py

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