<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Imports and setup
Imports
Step4:  Common functions
Step5:  Get charges
Calculate RESP charges using Gaussian through submit_gaussian for use with GAFF.
Step6:  Parameterize molecule in GAFF with ANTECHAMBER and ACPYPE
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.
For the catalyst
Step7:  Move molecules
In VMD, the molecules were moved so that they were not sitting on top of each other.
Solvate
As before, using DCM parameters and solvent box from virtualchemistry.org.
Step8:  Minimize
Step9:  Equilibrate
Step10:  Setup and submit parallel tempering (PT)
Step11:  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.
Step12:  Setup for several systems/molecules at once
Working based on what was done above (using some things that were defined up there as well
Get charges
Step13:  Copied over the g16.gesp files and renamed them for each molecule.
Make input files
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.
Step14:  Move molecules
I presume I will again need to make the molecules non-overlapping, and that will be done manually in VMD.
Box and solvate
Step15:  Minimize
Step16:  Made index file (called index-ycp.ndx) with solutes and solvent groups.
SA equilibration
Step17:  !!! Need to check distance on restraint !!!
Check equilibration
Step18:  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.
<ASSISTANT_TASK:>
Python Code:
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
Explanation: Imports and setup
Imports
End of explanation
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
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
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}')
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.')
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)
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)
Explanation: Common functions
End of explanation
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')
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
subg.submit_scripts(l_scripts, batch=True, submit=True)
Explanation: Get charges
Calculate RESP charges using Gaussian through submit_gaussian for use with GAFF.
End of explanation
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')
for res in gpx.residues:
    if res.name == 'MOL':
        res.name = 'GPX'
for res in cpa.residues:
    if res.name == 'MOL':
        res.name = 'CPA'
struc_comb = gpx + cpa
struc_comb
struc_comb.write('gpx-cpa-dry.top')
struc_comb.save('gpx-cpa-dry.gro')
Explanation: Parameterize molecule in GAFF with ANTECHAMBER and ACPYPE
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.
For the catalyst:
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
For the reactant:
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.
End of explanation
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))
Explanation: Move molecules
In VMD, the molecules were moved so that they were not sitting on top of each other.
Solvate
As before, using DCM parameters and solvent box from virtualchemistry.org.
End of explanation
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'
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)')
Explanation: Minimize
End of explanation
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
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)')
Explanation: Equilibrate
End of explanation
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,
                  )
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)
Explanation: Setup and submit parallel tempering (PT)
End of explanation
e_05s = dict()
for i in range(16):
    e_05s[i] = panedr.edr_to_df(f'05-PT/PT-out{i}.edr')
fig, axes = plot_prop_PT(e_05s, 'Pressure')
Explanation: 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.
End of explanation
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
subg.submit_scripts(l_scripts, batch=True, submit=True)
Explanation: Setup for several systems/molecules at once
Working based on what was done above (using some things that were defined up there as well
Get charges
End of explanation
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')
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'
ts1_en = ts1.copy(pmd.gromacs.GromacsTopologyFile)
ts3_en = ts3.copy(pmd.gromacs.GromacsTopologyFile)
ts1_en.coordinates = - ts1.coordinates
ts3_en.coordinates = - ts3.coordinates
sys_ts1 = ts1 + ycp
sys_ts1_en = ts1_en + ycp
sys_ts3 = ts3 + ycp
sys_ts3_en = ts3_en + ycp
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')
Explanation: Copied over the g16.gesp files and renamed them for each molecule.
Make input files
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.
End of explanation
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))
Explanation: Move molecules
I presume I will again need to make the molecules non-overlapping, and that will be done manually in VMD.
Box and solvate
End of explanation
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)
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)')
Explanation: Minimize
End of explanation
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
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)
Explanation: Made index file (called index-ycp.ndx) with solutes and solvent groups.
SA equilibration
End of explanation
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')
for key in l_syss:
    d = e_38s[key]
    fig, axes = plot_prop_PT(d, 'Volume')
Explanation: !!! Need to check distance on restraint !!!
Check equilibration
End of explanation
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]
Explanation: 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.
End of explanation
<END_TASK>