Source code for physical_validation.util.integrator

###########################################################################
#                                                                         #
#    physical_validation,                                                 #
#    a python package to test the physical validity of MD results         #
#                                                                         #
#    Written by Pascal T. Merz <pascal.merz@me.com>                       #
#               Michael R. Shirts <michael.shirts@colorado.edu>           #
#                                                                         #
#    Copyright (c) 2017-2021 University of Colorado Boulder               #
#              (c) 2012      The University of Virginia                   #
#                                                                         #
###########################################################################
r"""
This module contains low-level functionality of the
`physical_validation.integrator` module. The functions in this module should
generally not be called directly. Please use the high-level functions from
`physical_validation.integrator`.
"""
from typing import Callable, Dict, Optional, Tuple

import numpy as np

from . import plot


[docs]def calculate_rmsd( data: np.ndarray, time: Optional[np.ndarray] ) -> Tuple[float, float, float]: assert isinstance(data, np.ndarray) and data.ndim == 1 assert time is None or isinstance(time, np.ndarray) and time.ndim == 1 avg = float(data.mean()) if time is None: time = np.arange(data.size) fit = np.polyfit(time, data, 1) rmsd = float(data.std()) return avg, rmsd, float(fit[0])
[docs]def max_deviation(dts: np.ndarray, rmsds: np.ndarray) -> float: dt_ratio_2 = (dts[:-1] / dts[1:]) ** 2 rmsds = rmsds[:-1] / rmsds[1:] return np.max(np.abs(1 - rmsds / dt_ratio_2))
[docs]def check_convergence( const_traj: Dict[str, np.ndarray], convergence_test: Callable[[np.ndarray, np.ndarray], float], verbose: bool, screen: bool, filename: Optional[str], ) -> float: assert isinstance(const_traj, dict) assert len(const_traj) >= 2 if verbose: print("{:65s}".format("-" * 65)) print( "{:>10s} {:>10s} {:>10s} {:>10s} {:^21s}".format( "dt", "avg", "rmsd", "slope", "ratio" ) ) print("{:43s} {:>10s} {:>10s}".format("", "dt^2", "rmsd")) print("{:65s}".format("-" * 65)) prev = None results = {} for dt, traj in sorted(const_traj.items(), key=lambda x: float(x[0]), reverse=True): assert isinstance(traj, np.ndarray) assert traj.ndim == 1 or traj.ndim == 2 dt = float(dt) if traj.ndim == 1: data = traj time = None else: data = traj[1] time = traj[0] results[dt] = calculate_rmsd(data, time) if verbose: if prev is None: print( "{:10.4g} {:10.2f} {:10.2e} {:10.2e} {:>10s} {:>10s}".format( dt, results[dt][0], results[dt][1], results[dt][2], "--", "--" ) ) prev = [dt, results[dt][1]] else: print( "{:10.4g} {:10.2f} {:10.2e} {:10.2e} {:10.2f} {:10.2f}".format( dt, results[dt][0], results[dt][1], results[dt][2], prev[0] ** 2 / dt**2, prev[1] / results[dt][1], ) ) prev = [dt, results[dt][1]] if verbose: print("{:65s}".format("-" * 65)) dts = np.sort(np.array([float(dt) for dt in results.keys()]))[::-1] rmsds = np.array([float(results[dt][1]) for dt in dts]) do_plot = screen or filename is not None if do_plot: data = [ { "x": dts[1:], "y": rmsds[:-1] / rmsds[1:], "name": "Integrator convergence", "args": {"marker": "o"}, }, { "x": dts[1:], "y": (dts[:-1] / dts[1:]) ** 2, "name": "Expected convergence", }, ] max_y = max(np.max(data[0]["y"]), np.max(data[1]["y"])) plot.plot( data, legend="best", title="Actual vs. expected convergence", xlabel="Time step", ylabel="Convergence", xlim=(0, dts[1]), ylim=(0.5, max_y + 0.5), inv_x=True, filename=filename, screen=screen, ) return convergence_test(dts, rmsds)