In [ ]:
import re, os, sys, shutil
import shlex, subprocess
import glob
import pandas as pd
import panedr
import numpy as np
import MDAnalysis as mda
import nglview
import matplotlib.pyplot as plt
import parmed as pmd
import py
import scipy
from scipy import stats
from importlib import reload
from thtools import cd
from paratemp import copy_no_overwrite
from paratemp import geometries as gm
from paratemp import coordinate_analysis as ca
import paratemp.para_temp_setup as pts
import paratemp as pt
from gautools import submit_gaussian as subg
from gautools.tools import use_gen_template as ugt
In [ ]:
def plot_prop_PT(edict, prop):
fig, axes = plt.subplots(4, 4, figsize=(16,16))
for i in range(16):
ax = axes.flat[i]
edict[i][prop].plot(ax=ax)
fig.tight_layout()
return fig, axes
In [ ]:
def plot_e_props(df, labels, nrows=2, ncols=2):
fig, axes = plt.subplots(nrows, ncols, sharex=True)
for label, ax in zip(labels, axes.flat):
df[label].plot(ax=ax)
ax.set_title(label)
fig.tight_layout()
return fig, axes
In [ ]:
def plot_rd(univ): # rd = reaction distance
univ.calculate_distances(rd=(20,39))
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
univ.data.rd.plot(ax=axes[0])
univ.data.rd.hist(ax=axes[1], grid=False)
print(f'reaction distance mean: {univ.data.rd.mean():.2f} and sd: {univ.data.rd.std():.2f}')
return fig, axes
def plot_hist_dist(univ, name, indexes=None):
if indexes is not None:
kwargs = {name: indexes}
univ.calculate_distances(**kwargs)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
univ.data[name].plot(ax=axes[0])
univ.data[name].hist(ax=axes[1], grid=False)
print(f'{name} distance mean: {univ.data[name].mean():.2f} and sd: {univ.data[name].std():.2f}')
In [ ]:
def get_solvent_count_solvate(proc):
for line in proc.stdout.split('\n'):
m = re.search(r'(?:atoms\):\s+)(\d+)(?:\s+residues)', line)
if m:
return int(m.group(1))
else:
raise ValueError('Solvent count not found.')
In [ ]:
def set_solv_count(n_gro, s_count,
res_name='DCM', prepend='unequal-'):
"""
Remove solvent residues from the end of a gro file to match s_count
This assumes all non-solvent molecules are listed in the input gro
file before the solvent residues.
"""
bak_name = os.path.join(os.path.dirname(n_gro),
prepend+os.path.basename(n_gro))
copy_no_overwrite(n_gro, bak_name)
with open(n_gro, 'r') as in_gro:
lines = in_gro.readlines()
for line in lines[2:]:
if res_name in line:
non_s_res_count = resid
break
else:
resid = int(line[:5])
res_count = s_count + non_s_res_count
# TODO check reasonability of this number
box = lines.pop()
while True:
line = lines.pop()
if int(line[:5]) > res_count:
continue
elif int(line[:5]) == res_count:
atom_count = line[15:20]
lines.append(line)
break
elif int(line[:5]) < res_count:
raise ValueError("Desired res "
"count is larger than "
"line's resid.\n" +
"res_count: {}\n".format(res_count) +
"line: {}".format(line))
lines[1] = atom_count + '\n'
lines.append(box)
with open(n_gro, 'w') as out_gro:
for line in lines:
out_gro.write(line)
In [ ]:
def get_solv_count_top(n_top, res_name='DCM'):
"""
Return residue count of specified residue from n_top"""
with open(n_top, 'r') as in_top:
mol_section = False
for line in in_top:
if line.strip().startswith(';'):
pass
elif not mol_section:
if re.search(r'\[\s*molecules\s*\]', line,
flags=re.IGNORECASE):
mol_section = True
else:
if res_name.lower() in line.lower():
return int(line.split()[1])
def set_solv_count_top(n_top, s_count,
res_name='DCM', prepend='unequal-'):
"""
Set count of res_name residues in n_top
This will make a backup copy of the top file with `prepend`
prepended to the name of the file."""
bak_name = os.path.join(os.path.dirname(n_top),
prepend+os.path.basename(n_top))
copy_no_overwrite(n_top, bak_name)
with open(n_top, 'r') as in_top:
lines = in_top.readlines()
with open(n_top, 'w') as out_top:
mol_section = False
for line in lines:
if line.strip().startswith(';'):
pass
elif not mol_section:
if re.search(r'\[\s*molecules\s*\]', line,
flags=re.IGNORECASE):
mol_section = True
else:
if res_name.lower() in line.lower():
line = re.sub(r'\d+', str(s_count), line)
out_top.write(line)
Calculate RESP charges using Gaussian through submit_gaussian for use with GAFF.
In [ ]:
d_charge_params = dict(opt='SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) iop(6/50=1)',
func='HF',
basis='6-31G*',
footer='\ng16.gesp\n\ng16.gesp\n\n')
In [ ]:
l_scripts = []
s = subg.write_sub_script('01-charges/TS2.com',
executable='g16',
make_xyz='../TS2.pdb',
make_input=True,
ugt_dict={'job_name':'GPX TS2 charges',
'charg_mult':'+1 1',
**d_charge_params})
l_scripts.append(s)
s = subg.write_sub_script('01-charges/R-NO2-CPA.com',
executable='g16',
make_xyz='../R-NO2-CPA.pdb',
make_input=True,
ugt_dict={'job_name':'GPX R-NO2-CPA charges',
'charg_mult':'-1 1',
**d_charge_params})
l_scripts.append(s)
l_scripts
In [ ]:
subg.submit_scripts(l_scripts, batch=True, submit=True)
Note, ACPYPE was installed from this repository, which seems to be from the original author, though maybe not the one who put it onto pypi.
Use antechamber to create mol2 file with Gaussian ESP charges (though wrong atom types and such, for now):
antechamber -i R-NO2-CPA.gesp -fi gesp -o R-NO2-CPA.mol2 -fo mol2
Use ACPYPE to use this mol2 file (and it's GESP charges) to generate GROMACS input files:
acpype.py -i R-NO2-CPA.mol2 -b CPA-gesp --net_charge=-1 -o gmx -d -c user
antechamber -i TS2.gesp -fi gesp -o TS2.mol2 -fo mol2
acpype.py -i TS2.mol2 -b GPX-ts --net_charge=1 -o gmx -c user
Then the different molecules can be combined using ParmEd.
In [ ]:
gpx = pmd.gromacs.GromacsTopologyFile('01-charges/GPX-ts.acpype/GPX-ts_GMX.top', xyz='01-charges/GPX-ts.acpype/GPX-ts_GMX.gro')
cpa = pmd.gromacs.GromacsTopologyFile('01-charges/CPA-gesp.acpype/CPA-gesp_GMX.top', xyz='01-charges/CPA-gesp.acpype/CPA-gesp_GMX.gro')
In [ ]:
for res in gpx.residues:
if res.name == 'MOL':
res.name = 'GPX'
for res in cpa.residues:
if res.name == 'MOL':
res.name = 'CPA'
In [ ]:
struc_comb = gpx + cpa
struc_comb
In [ ]:
struc_comb.write('gpx-cpa-dry.top')
In [ ]:
struc_comb.save('gpx-cpa-dry.gro')
In [ ]:
f_dcm = py.path.local('~/GROMACS-basics/DCM-GAFF/')
f_solvate = py.path.local('02-solvate/')
sep_gro = py.path.local('gpx-cpa-sep.gro')
boxed_gro = f_solvate.join('gpx-cpa-boxed.gro')
box = '3.5 3.5 3.5'
solvent_source = f_dcm.join('dichloromethane-T293.15.gro')
solvent_top = f_dcm.join('dichloromethane.top')
solv_gro = f_solvate.join('gpx-cpa-dcm.gro')
top = py.path.local('../params/gpxTS-cpa-dcm.top')
verbose = True
solvent_counts, key = dict(), 'GPX'
with f_solvate.as_cwd():
## Make box
cl = shlex.split(f'gmx_mpi editconf -f {sep_gro} ' +
f'-o {boxed_gro} -box {box}')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_editconf'] = proc.stdout
proc.check_returncode()
## Solvate
cl = shlex.split(f'gmx_mpi solvate -cp {boxed_gro} ' +
f'-cs {solvent_source} -o {solv_gro}')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_solvate'] = proc.stdout
proc.check_returncode()
solvent_counts[key] = get_solvent_count_solvate(proc)
if verbose:
print(f'Solvated system into {solv_gro}')
struc_g_c = pmd.load_file('gpx-cpa-dry.top')
struc_dcm = pmd.load_file(str(f_dcm.join('dichloromethane.top')))
struc_g_c_d = struc_g_c + solvent_counts['GPX'] * struc_dcm
struc_g_c_d.save(str(top))
In [ ]:
ppl = py.path.local
f_min = ppl('03-minimize/')
f_g_basics = py.path.local('~/GROMACS-basics/')
mdp_min = f_g_basics.join('minim.mdp')
tpr_min = f_min.join('min.tpr')
deffnm_min = f_min.join('min-out')
gro_min = deffnm_min + '.gro'
In [ ]:
with f_min.as_cwd():
## Compile tpr
if not tpr_min.exists():
cl = shlex.split(f'gmx_mpi grompp -f {mdp_min} '
f'-c {solv_gro} '
f'-p {top} '
f'-o {tpr_min}')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_grompp_em'] = proc.stdout
proc.check_returncode()
if verbose:
print(f'Compiled em tpr to {tpr_min}')
elif verbose:
print(f'em tpr file already exists ({tpr_min})')
## Run minimization
if not gro_min.exists():
cl = shlex.split('gmx_mpi mdrun '
f'-s {tpr_min} '
f'-deffnm {deffnm_min} ')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_mdrun_em'] = proc.stdout
# TODO Get the potential energy from this output
proc.check_returncode()
if verbose:
print(f'Ran {key} em to make {gro_min}')
elif verbose:
print(f'em output gro already exists (gro_min)')
In [ ]:
f_equil = ppl('04-equilibrate/')
plumed = f_equil.join('plumed.dat')
mdp_equil = f_g_basics.join('npt-298.mdp')
tpr_equil = f_equil.join('equil.tpr')
deffnm_equil = f_equil.join('equil-out')
gro_equil = deffnm_equil + '.gro'
gro_input = gro_min
In [ ]:
with f_equil.as_cwd():
## Compile equilibration
if not tpr_equil.exists():
cl = shlex.split(f'gmx_mpi grompp -f {mdp_equil} '
f'-c {gro_input} '
f'-p {top} '
f'-o {tpr_equil}')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_grompp_equil'] = proc.stdout
proc.check_returncode()
if verbose:
print(f'Compiled equil tpr to {tpr_equil}')
elif verbose:
print(f'equil tpr file already exists ({tpr_equil})')
## Run equilibration
if not gro_equil.exists():
cl = shlex.split('gmx_mpi mdrun '
f'-s {tpr_equil} '
f'-deffnm {deffnm_equil} '
f'-plumed {plumed}')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_mdrun_equil'] = proc.stdout
proc.check_returncode()
if verbose:
print(f'Ran {key} equil to make {gro_equil}')
elif verbose:
print(f'equil output gro already exists (gro_equil)')
In [ ]:
f_pt = ppl('05-PT/')
template = f_pt.join('template-mdp.txt')
index = ppl('index.ndx')
sub_templ = f_g_basics.join('sub-template-128.sub')
d_sub_templ = dict(tpr_base = 'TOPO/npt',
deffnm = 'PT-out',
name = 'GPX-PT',
plumed = plumed,
)
In [ ]:
scaling_exponent = 0.025
maxwarn = 0
start_temp = 298.
verbose = True
skip_existing = True
jobs = []
failed_procs = []
for key in ['GPX']:
kwargs = {'template': str(template),
'topology': str(top),
'structure': str(gro_equil),
'index': str(index),
'scaling_exponent': scaling_exponent,
'start_temp': start_temp,
'maxwarn': maxwarn}
with f_pt.as_cwd():
try:
os.mkdir('TOPO')
except FileExistsError:
if skip_existing:
print(f'Skipping {key} because it seems to '
'already be done.\nMoving on...')
continue
with cd('TOPO'):
print(f'Now in {os.getcwd()}\nAttempting to compile TPRs...')
pts.compile_tprs(**kwargs)
print('Done compiling. Moving on...')
print(f'Now in {os.getcwd()}\nWriting submission script...')
with sub_templ.open(mode='r') as templ_f, \
open('gromacs-start-job.sub', 'w') as sub_s:
[sub_s.write(l.format(**d_sub_templ)) for l in templ_f]
print('Done.\nNow submitting job...')
cl = ['qsub', 'gromacs-start-job.sub']
proc = subprocess.run(cl,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT,
universal_newlines=True)
if proc.returncode == 0:
output = proc.stdout
jobs.append(re.search('[0-9].+\)', output).group(0))
print(output, '\nDone.\nMoving to next...')
else:
print('\n\n'+5*'!!!---'+'\n')
print(f'Error with calling qsub on {key}')
print('Command line input was', cl)
print('Check input and try again manually.'
'\nMoving to next anyway...')
failed_procs.append(proc)
print('-----Done-----\nSummary of jobs submitted:')
for job in jobs:
print(job)
The energies from the simulations can be read in as a pandas DataFrame using panedr and then analyzed or plotted to check on equilibration, convergence, etc.
In [ ]:
e_05s = dict()
for i in range(16):
e_05s[i] = panedr.edr_to_df(f'05-PT/PT-out{i}.edr')
In [ ]:
fig, axes = plot_prop_PT(e_05s, 'Pressure')
In [ ]:
l_scripts = []
s = subg.write_sub_script('01-charges/TS1.com',
executable='g16',
make_xyz='../TS1protonated.mol2',
make_input=True,
ugt_dict={'job_name':'GPX TS1 charges',
'charg_mult':'+1 1',
**d_charge_params})
l_scripts.append(s)
s = subg.write_sub_script('01-charges/TS3.com',
executable='g16',
make_xyz='../TS3protonated.mol2',
make_input=True,
ugt_dict={'job_name':'GPX TS3 charges',
'charg_mult':'+1 1',
**d_charge_params})
l_scripts.append(s)
s = subg.write_sub_script('01-charges/anti-cat-yamamoto.com',
executable='g16',
make_xyz='../R-Yamamoto-Cat.pdb',
make_input=True,
ugt_dict={'job_name':
'yamamoto catalyst charges',
'charg_mult':'-1 1',
**d_charge_params})
l_scripts.append(s)
l_scripts
In [ ]:
subg.submit_scripts(l_scripts, batch=True, submit=True)
Copied over the g16.gesp files and renamed them for each molecule.
Loaded amber/2016 module (and its dependencies).
antechamber -i TS1.gesp -fi gesp -o TS1.mol2 -fo mol2
acpype.py -i TS1.mol2 -b TS1-gesp --net_charge=1 -o gmx -d -c user
There was a warning for assigning bond types.
antechamber -i TS3.gesp -fi gesp -o TS3.mol2 -fo mol2
acpype.py -i TS3.mol2 -b TS3-gesp --net_charge=1 -o gmx -d -c user
Similar warning.
antechamber -i YCP.gesp -fi gesp -o YCP.mol2 -fo mol2
acpype.py -i YCP.mol2 -b YCP-gesp --net_charge=-1 -o gmx -d -c use
No similar warning here.
In [ ]:
ts1 = pmd.gromacs.GromacsTopologyFile(
'01-charges/TS1-gesp.acpype/TS1-gesp_GMX.top',
xyz='01-charges/TS1-gesp.acpype/TS1-gesp_GMX.gro')
ts3 = pmd.gromacs.GromacsTopologyFile(
'01-charges/TS3-gesp.acpype/TS3-gesp_GMX.top',
xyz='01-charges/TS3-gesp.acpype/TS3-gesp_GMX.gro')
ycp = pmd.gromacs.GromacsTopologyFile(
'01-charges/YCP-gesp.acpype/YCP-gesp_GMX.top',
xyz='01-charges/YCP-gesp.acpype/YCP-gesp_GMX.gro')
In [ ]:
for res in ts1.residues:
if res.name == 'MOL':
res.name = 'TS1'
for res in ts3.residues:
if res.name == 'MOL':
res.name = 'TS3'
for res in ycp.residues:
if res.name == 'MOL':
res.name = 'YCP'
In [ ]:
ts1_en = ts1.copy(pmd.gromacs.GromacsTopologyFile)
ts3_en = ts3.copy(pmd.gromacs.GromacsTopologyFile)
ts1_en.coordinates = - ts1.coordinates
ts3_en.coordinates = - ts3.coordinates
In [ ]:
sys_ts1 = ts1 + ycp
sys_ts1_en = ts1_en + ycp
sys_ts3 = ts3 + ycp
sys_ts3_en = ts3_en + ycp
In [ ]:
sys_ts1.write('ts1-ycp-dry.top')
sys_ts3.write('ts3-ycp-dry.top')
sys_ts1.save('ts1-ycp-dry.gro')
sys_ts1_en.save('ts1_en-ycp-dry.gro')
sys_ts3.save('ts3-ycp-dry.gro')
sys_ts3_en.save('ts3_en-ycp-dry.gro')
In [ ]:
f_dcm = py.path.local('~/GROMACS-basics/DCM-GAFF/')
f_solvate = py.path.local('37-solvate-anti/')
box = '3.7 3.7 3.7'
solvent_source = f_dcm.join('dichloromethane-T293.15.gro')
solvent_top = f_dcm.join('dichloromethane.top')
solv_gro = f_solvate.join('gpx-cpa-dcm.gro')
ts1_top = ppl('../params/ts1-ycp-dcm.top')
ts3_top = ppl('../params/ts3-ycp-dcm.top')
l_syss = ['TS1', 'TS1_en', 'TS3', 'TS3_en']
verbose = True
solvent_counts = dict()
for key in l_syss:
sep_gro = ppl(f'{key.lower()}-ycp-dry.gro')
if not sep_gro.exists():
raise FileNotFoundError(f'{sep_gro} does not exist')
boxed_gro = f'{key.lower()}-ycp-box.gro'
solv_gro = f'{key.lower()}-ycp-dcm.gro'
with f_solvate.ensure_dir().as_cwd():
## Make box
cl = shlex.split(f'gmx_mpi editconf -f {sep_gro} ' +
f'-o {boxed_gro} -box {box}')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_editconf'] = proc.stdout
proc.check_returncode()
## Solvate
cl = shlex.split(f'gmx_mpi solvate -cp {boxed_gro} ' +
f'-cs {solvent_source} -o {solv_gro}')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_solvate'] = proc.stdout
proc.check_returncode()
solvent_counts[key] = get_solvent_count_solvate(proc)
if verbose:
print(f'Solvated system into {solv_gro}')
# min_solv_count = min(solvent_counts.values())
min_solv_count = 328 # want to match with syn calculations
if min(solvent_counts.values()) < min_solv_count:
raise ValueError('At least one of the structures has <328 DCMs.\n'
'Check and/or make the box larger')
for key in l_syss:
solv_gro = f'{key.lower()}-ycp-dcm.gro'
with f_solvate.as_cwd():
set_solv_count(solv_gro, min_solv_count)
struc_ts1 = pmd.load_file('ts1-ycp-dry.top')
struc_ts3 = pmd.load_file('ts3-ycp-dry.top')
struc_dcm = pmd.load_file(str(f_dcm.join('dichloromethane.top')))
struc_ts1_d = struc_ts1 + min_solv_count * struc_dcm
struc_ts1_d.save(str(ts1_top))
struc_ts3_d = struc_ts3 + min_solv_count * struc_dcm
struc_ts3_d.save(str(ts3_top))
In [ ]:
f_min = ppl('38-relax-anti/')
f_min.ensure_dir()
f_g_basics = py.path.local('~/GROMACS-basics/')
mdp_min = f_g_basics.join('minim.mdp')
d_tops = dict(TS1=ts1_top, TS1_en=ts1_top, TS3=ts3_top, TS3_en=ts3_top)
In [ ]:
for key in l_syss:
solv_gro = ppl(f'37-solvate-anti/{key.lower()}-ycp-dcm.gro')
tpr_min = f_min.join(f'{key.lower()}-min.tpr')
deffnm_min = f_min.join(f'{key.lower()}-min-out')
gro_min = deffnm_min + '.gro'
top = d_tops[key]
with f_min.as_cwd():
## Compile tpr
if not tpr_min.exists():
cl = shlex.split(f'gmx_mpi grompp -f {mdp_min} '
f'-c {solv_gro} '
f'-p {top} '
f'-o {tpr_min}')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_grompp_em'] = proc.stdout
proc.check_returncode()
if verbose:
print(f'Compiled em tpr to {tpr_min}')
elif verbose:
print(f'em tpr file already exists ({tpr_min})')
## Run minimization
if not gro_min.exists():
cl = shlex.split('gmx_mpi mdrun '
f'-s {tpr_min} '
f'-deffnm {deffnm_min} ')
proc = subprocess.run(cl, universal_newlines=True,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT)
outputs[key+'_mdrun_em'] = proc.stdout
# TODO Get the potential energy from this output
proc.check_returncode()
if verbose:
print(f'Ran {key} em to make {gro_min}')
elif verbose:
print(f'em output gro already exists (gro_min)')
In [ ]:
f_pt = ppl('38-relax-anti/')
template = ppl('33-SA-NPT-rest-no-LINCS/template-mdp.txt')
index = ppl('../params/index-ycp.ndx')
scaling_exponent = 0.025
maxwarn = 0
start_temp = 298.
nsims = 16
In [ ]:
verbose = True
skip_existing = True
jobs = []
failed_procs = []
for key in l_syss:
d_sub_templ = dict(
tpr = f'{key.lower()}-TOPO/npt',
deffnm = f'{key.lower()}-SA-out',
name = f'{key.lower()}-SA',
nsims = nsims,
tpn = 16,
cores = 128,
multi = True,
)
gro_equil = f_min.join(f'{key.lower()}-min-out.gro')
top = d_tops[key]
kwargs = {'template': str(template),
'topology': str(top),
'structure': str(gro_equil),
'index': str(index),
'scaling_exponent': scaling_exponent,
'start_temp': start_temp,
'maxwarn': maxwarn,
'number': nsims,
'grompp_exe': 'gmx_mpi grompp'}
with f_pt.as_cwd():
try:
os.mkdir(f'{key.lower()}-TOPO/')
except FileExistsError:
if (os.path.exists(f'{key.lower()}-TOPO/temperatures.dat') and
skip_existing):
print(f'Skipping {key} because it seems to '
'already be done.\nMoving on...')
continue
with cd(f'{key.lower()}-TOPO/'):
print(f'Now in {os.getcwd()}\nAttempting to compile TPRs...')
pts.compile_tprs(**kwargs)
print('Done compiling. Moving on...')
print(f'Now in {os.getcwd()}\nWriting submission script...')
lp_sub = pt.sim_setup.make_gromacs_sub_script(
f'gromacs-start-{key}-job.sub', **d_sub_templ)
print('Done.\nNow submitting job...')
cl = shlex.split(f'qsub {lp_sub}')
proc = subprocess.run(cl,
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT,
universal_newlines=True)
if proc.returncode == 0:
output = proc.stdout
jobs.append(re.search('[0-9].+\)', output).group(0))
print(output, '\nDone.\nMoving to next...')
else:
print('\n\n'+5*'!!!---'+'\n')
print(f'Error with calling qsub on {key}')
print('Command line input was', cl)
print('Check input and try again manually.'
'\nMoving to next anyway...')
failed_procs.append(proc)
print('-----Done-----\nSummary of jobs submitted:')
for job in jobs:
print(job)
!!! Need to check distance on restraint !!!
In [ ]:
e_38s = dict()
for key in l_syss:
deffnm = f'{key.lower()}-SA-out'
e_38s[key] = dict()
d = e_38s[key]
for i in range(16):
d[i] = panedr.edr_to_df(f'38-relax-anti/{deffnm}{i}.edr')
In [ ]:
for key in l_syss:
d = e_38s[key]
fig, axes = plot_prop_PT(d, 'Volume')
The volumes seem to look okay. Started high (I did remove some solvents and it hadn't relaxed much), dropped quickly, then seemed to grow appropriately as the temperatures rose. None seems to have boiled.
In [ ]:
for key in l_syss:
d = e_38s[key]
fig, ax = plt.subplots()
for key in list(d.keys()):
ax.hist(d[key]['Total Energy'], bins=100)
del d[key]
In [ ]: