Source code for schrodinger.application.desmond.deltaE

import os
import re
from collections import defaultdict

import numpy as np


def gen_dE_files(fname, out_fname):
    """
    'fname' is the output deltaE file.
    'out_fname' is a dict: Key = replica index (0-based); value = gibbs' .dE file name.
    return True if successful; return False if any output file is empty
    """
    REPLICA = re.compile(r"\[T(\d+),T(\d+)\] *: *\[(.+),(.+)\]")
    # Key = (host-replica-index, trial-replica-index); value = List of floating numbers (energy)
    energy = defaultdict(list)
    timeline = []
    with open(fname, "r") as fh:
        for line in fh.readlines():
            line = line.strip()
            if 'Time' == line[:4]:
                timeline.append(float(line[7:]))
            elif '[' == line[0]:
                match = REPLICA.match(line)
                if match:
                    r0, r1, ef, er = match.group(1, 2, 3, 4)
                    r0, r1 = map(int, (r0, r1))
                    ef, er = map(float, (ef, er))
                    energy[(r0, r1)].append(ef)
                    energy[(r1, r0)].append(er)

    # it is possible the files are truncated, for instance when monitoring
    # convergence of running simulation
    minsize = min([len(e) for e in energy.values()])

    all_found = True
    for i in range(len(out_fname)):
        ef = energy[(i, i + 1)][:minsize]
        er = energy[(i, i - 1)][:minsize]
        timeline = timeline[:minsize]

        if not (er or ef):
            all_found = False
            continue

        ef = ef or np.zeros_like(er)
        er = er or np.zeros_like(ef)

        dE_array = np.stack([timeline, er, ef], axis=-1)
        fname_base = os.path.basename(fname)
        additional_header = f"Output recreated from REMD {fname_base} file"
        write_dE(out_fname[i], dE_array, additional_header=additional_header)

    return all_found


def write_dE(out_fname: str, dE_array: np.array, additional_header: str = ""):
    """
    :param out_fname: The output filename
    :param dE_array: An array with shape (Nx3) where N is the number of output
        timesteps with columns representing the time, reverse dE, and forward
        dE.
    :param additional_header: Header text which will be prepended before the
        column headerss
    """
    header = "\n".join(
        (additional_header, "0:time (ps)  1:dE- (kcal/mol)  2:dE+ (kcal/mol)"))
    format = ("%13.3f", "%18.8e", "%18.8e")
    np.savetxt(out_fname, dE_array, header=header, fmt=format, comments='# ')