import os
import shlex
import subprocess
import tempfile
from past.utils import old_div
import numpy
[docs]def generate_config(ensemble,
bootfile,
lastime=10.,
gpu=True,
dsmworkdir=None,
additional_config=None):
# Generate either a configuration file that can be passed to desmond(gpu),
# or create a run directory that can be passed to dsm_[sbatch,salloc]
options = []
# Select ensemble
if ensemble != "NPT":
options += '''
--cfg integrator.Multigrator.barostat=none
'''.split()
if ensemble != "NVT":
options += '''
--cfg integrator.Multigrator.thermostat=none
'''.split()
# Switch to Verlet
options += shlex.split('''
--cfg integrator.Multigrator.nve.type=Verlet
--cfg integrator.respa.far_timesteps=1
--cfg integrator.dt=0.0025
--cfg "force.nonbonded.far.order=[ 6 6 6 ]"
''')
# Set simulation time
options += shlex.split('''
--cfg mdsim.last_time=%g
''' % lastime)
# Default output
options += shlex.split('''
--cfg "mdsim.plugin.list[+]=eneseq"
--cfg "mdsim.plugin.list[+]=status"
''')
if gpu:
options += ['--gpu']
retval = None
if dsmworkdir is not None:
# Set up for use with dms_sbatch
retval = dsmworkdir
options += shlex.split('''
--dms-configure %s''' % retval)
# Do not create a locker for our tests
options += shlex.split('''
--cfg desmond_submit.locker.enable=false
--cfg desmond_submit.modules=[]
''')
else:
# Set up for local run
retval = tempfile.mktemp()
options += shlex.split('''
--output-ark %s''' % retval)
if additional_config is not None:
options += additional_config
command = shlex.split('''
desmond_configure_for.py %s
''' % (bootfile))
command += options
subprocess.check_output(command, universal_newlines=True)
return retval
[docs]def parse_eneseq(filename):
# Read an eneseq file and return its data as a dict which
# maps from column name to column data
column_to_key = [
'time', 'conserved energy', 'potential energy', 'kinetic energy',
'com energy', 'extended energy', 'force energy', 'pressure', 'volume',
'temperature'
]
edict = dict()
data = numpy.genfromtxt(filename)
if len(data.shape) == 1:
data = data.reshape([1, data.size])
for i, k in enumerate(column_to_key):
edict[k] = data[:, i]
return edict
[docs]def compare_eneseq(edict1, edict2):
# Compare two parsed eneseq structure and report the
# mismatches based on the hard-coded tolerances.
tolerances = {
'time': None,
'conserved energy': (1.e-3, 1.e-5),
'potential energy': (1.e-3, 1.e-5),
'kinetic energy': (1.e-3, 1.e-5),
'com energy': None,
'extended energy': None,
'force energy': (1.e-3, 1.e-5),
'pressure': (1.0e-1, 1.0e-4),
'volume': None,
'temperature': (1.0e-2, 1.0e-4)
}
times1 = edict1["time"]
times2 = edict2["time"]
if times1.size > times2.size:
times = times1
else:
times = times2
failed = []
for key, tol in tolerances.items():
if tol is None:
continue
ecol1 = edict1[key]
ecol2 = edict2[key]
diffa = (ecol1 - ecol2)
diffr = old_div(diffa, (ecol1 + 1.0e-30))
conda = abs(diffa) > tol[0]
condr = abs(diffr) > tol[1]
error = numpy.append([diffa], [diffr], 0).min(0)
failures = numpy.nonzero(numpy.all([conda, condr], axis=0))[0]
if len(failures) > 0:
failed.append(key)
print(key)
print("-" * 80)
for f in failures:
print(
"time: %12.6f, 1: %20.12e 2: %20.12e, errors: %g (a), %g (r)"
% (times[f], ecol1[f], ecol2[f], diffa[f], diffr[f]))
print(" ")
return failed
[docs]def run_desmond(configfile, gpu=True, additional_config=None, verbose=False):
eneseqfile = tempfile.mktemp() + ".eneseq"
executable = None
if gpu:
executable = os.path.join(os.environ["GDESMOND_PREFIX"], "bin/gdesmond")
else:
executable = os.path.join(os.environ["DESMOND_PREFIX"], "bin/desmond")
command = [executable]
if not gpu:
command += ['--tpp', '8']
command += shlex.split(" --include %s" % configfile)
command += shlex.split(" --cfg mdsim.plugin.eneseq.name=%s" % eneseqfile)
if additional_config is not None:
command += additional_config
if verbose:
print("# ", " ".join(command))
subprocess.check_output(command,
stderr=open("/tmp/gdesmond.stderr", "w"),
universal_newlines=True)
return eneseqfile