Strychnine

Along the steps of this example workflow we will show how to:

  1. Generate different conformers of the molecule using CSEARCH

  2. Generate the inputs for the QM geometry optimization

  3. Fix error terminations and imaginary frequencies of the output files

  4. Calculation and analyze the NMR chemical shifts for the conformers generated.

  5. Use GoodVibes to calculate the Boltzmann distributions using Gibbs free energies at 298.15 K

Specifically, in this example we will calculate the NMR chemical shifts of the strychnine starting from the smiles representation of said molecule that we can see below.

SMILES

C1CN2CC3=CCO[C@H]4CC(=O)N5[C@H]6[C@H]4[C@H]3C[C@H]2[C@@]61C7=CC=CC=C75

Strychnine

Note

A jupyter notebook containing all the steps shown in this example can be found in the AQME repository in Github or in Figshare

Note

A video tutorial illustrating this example can be found in our youtube channel

Step 1: Importing AQME and other python modules

import os, subprocess,shutil, glob
from pathlib import Path
from aqme.csearch import csearch
from aqme.qprep import qprep
from aqme.qcorr import qcorr
from aqme.qdescp import qdescp

Step 2: CSEARCH conformational sampling

name = 'Strychnine'
smi = 'C1CN2CC3=CCO[C@H]4CC(=O)N5[C@H]6[C@H]4[C@H]3C[C@H]2[C@@]61C7=CC=CC=C75'
program = 'rdkit'

# folder where the SDF files are generated
sdf_path = f'{os.getcwd()}/{name}_sdf-files'

csearch(destination=sdf_path,
        program=program,
        smi=smi,
        name=name)

Step 3: Creating Gaussian input files for optimization and frequency with QPREP

program = 'gaussian'
qm_input = 'B3LYP/6-31+G(d,p) opt freq'
mem='24GB'
nprocs=12

# SDF files from Step 2
sdf_rdkit_files = f'{sdf_path}/*.sdf'

# folder where the COM files are generated
com_path = f'{os.getcwd()}/{name}_com-files'

qprep(destination=com_path,
      files=sdf_rdkit_files,
      program=program,
      qm_input=qm_input,
      mem=mem,
      nprocs=nprocs)

Step 4: Running Gaussian inputs for optimization and frequency calcs externally

Now that we have generated our gaussian input files (in the com_path location of Step 3) we need to run the gaussian calculations. If you do not know how to run the Gaussian calculations in your HPC please refer to your HPC manager.

As an example, for a single calculation in Gaussian 16 through the terminal we would run the following command on a Linux-based system:

g16 myfile.com

Step 5: QCORR analysis including isomerization filter

log_files=f'{com_path}/*.log' # LOG files from Step 4

qcorr(files=log_files,
      freq_conv='opt=(calcfc,maxstep=5)',
      isom_type='com',
      isom_inputs=com_path,
      nprocs=12,
      mem='24GB')

Step 6: Resubmission of unsuccessful calculations (if any) with suggestions from AQME

Now we need to run the generated COM files (in fixed_QM_inputs) with Gaussian like we did in Step 4

After the calculations finish we check again the files using QCORR

new_log_files = f'{com_path}/failed/run_1/fixed_QM_inputs/*.log'

qcorr(files=new_log_files,
       isom_type='com',
       isom_inputs=f'{com_path}/failed/run_1/fixed_QM_inputs',
       nprocs=12,
       mem='24GB')

Step 7: Creating Gaussian input files for NMR calcs with QPREP

program = 'gaussian'
qm_input = 'B3LYP/6-311+G(2d,p) scrf=(solvent=chloroform,smd) nmr=giao'
mem='24GB'
nprocs=12

# folder where the successful LOG files are stored during the QCORR cycles
# (Steps 5 and 6)

success_folder = com_path+'/success'
log_files = f'{success_folder}/*.log'

# folder to store the new COM inputs for single point NMR calcs
sp_path = f'{os.getcwd()}/{name}_sp-files'

qprep(w_dir_main=success_folder,
      destination=sp_path,
      files=log_files,
      program=program,
      qm_input=qm_input,
      mem=mem,
      nprocs=nprocs,
      suffix='SP')

Step 8: Running Gaussian NMR calcs

Now we need to run the generated COM files (in sp_path) with Gaussian like we did in Step 4

After the calculations end, we create JSON files with QCORR to store the information from the resulting LOG files

log_files=f'{sp_path}/*.log'
qcorr(files=log_files)

Step 9: Obtaining Boltzmann weighted NMR shifts with QDESCP

# Analyze the JSON files to calculate the Boltzmann averaged shielding tensors

## folder where the JSON files were just created with QCORR
json_folder = sp_path+'/success/SP_calcs/json_files'
json_files=f'{json_folder}/*.json'

## folder to store the results from QDESCP
nmr_path = f'{os.getcwd()}/{name}_nmr-files'

qdescp(program='nmr',
       boltz=True,
       files=json_files,
       destination=nmr_path,
       nmr_slope=[-1.0537, -1.0784],
       nmr_intercept=[181.7815,31.8723],
       nmr_experim='Experimental_NMR_shifts.csv')

Step 10: Calculating conformer populations with GoodVibes

opt_files = glob.glob(f'{success_folder}/*.log')
sp_files = glob.glob(f'{sp_path}/success/SP_calcs/*.log')
log_files = opt_files + sp_files

w_dir_main  = Path(os.getcwd())
GV_folder = w_dir_main.joinpath('Strychine_GoodVibes-analysis')
GV_folder.mkdir(exist_ok=True, parents=True)

for file in log_files:
    shutil.copy(file, GV_folder)

# run GoodVibes
os.chdir(GV_folder)
subprocess.run(['python', '-m', 'goodvibes', '--xyz','-c','1', '*.log','--boltz'])
os.chdir(w_dir_main)