Quickstart guide with example input file

A short example illustrating how to run the complete workflow, and how to restart from previous calculations.

#!/bin/sh

# Reactions discovery example for a mixture of CH3CH2 and OH radicals
# The mixture should form ethanol CH3CH2OH and hydrogen peroxide HOOH
# and possibly more products.

rm -rf complete_nanoreactor.results

AMS_JOBNAME=complete_nanoreactor $AMSBIN/reactions_discovery << eor
   MolecularDynamics
      Enabled Yes
      NumSimulations 4
      BuildSystem
         NumAtoms 50
         Density 0.9
         Molecule 
            SMILES C[CH2]  # a CH3CH2 radical
            MoleFraction 1
         End
         Molecule 
            SMILES [OH]   # an OH radical
            MoleFraction 3
         End
      End
      Type NanoReactor
      NanoReactor
         NumCycles 5
         Temperature 500
         MinVolumeFraction 0.6
      End
      BondOrders Method=Guess
   End

   NetworkExtraction Enabled=Yes UseCharges=Yes
   ProductRanking Enabled=Yes

   Engine ReaxFF
      ForceField CHON-2019.ff
      TaperBO Yes
   EndEngine
eor

#Test if ethanol is found:
echo ETHANOL: `grep -c " CCO " complete_nanoreactor.results/reactions_discovery.log`
#Test if hydrogen peroxide is found:
echo HYDROGEN PEROXIDE: `grep -c " OO " complete_nanoreactor.results/reactions_discovery.log`

# Use the MD trajectories from the previous job
rm -rf restart_extraction.results
AMS_JOBNAME=restart_extraction $AMSBIN/reactions_discovery << eor
    MolecularDynamics Enabled=No
    NetworkExtraction Enabled=Yes MDTrajectories=complete_nanoreactor.results
    ProductRanking Enabled=Yes
    Engine ReaxFF
        ForceField CHON-2019.ff
        TaperBO Yes
    EndEngine
eor

#Test if ethanol is found:
echo ETHANOL: `grep -c " CCO " restart_extraction.results/reactions_discovery.log`
#Test if hydrogen peroxide is found:
echo HYDROGEN PEROXIDE: `grep -c " OO " restart_extraction.results/reactions_discovery.log`

# Use the previous NetworkExtraction results to restart the ranking
rm -rf restart_ranking.results
AMS_JOBNAME=restart_ranking $AMSBIN/reactions_discovery << eor
    MolecularDynamics Enabled=No
    NetworkExtraction Enabled=No
    ProductRanking Enabled=Yes ReactionNetwork=restart_extraction.results
eor

#Test if ethanol is found:
echo ETHANOL: `grep -c " CCO " restart_ranking.results/reactions_discovery.log`
#Test if hydrogen peroxide is found:
echo HYDROGEN PEROXIDE: `grep -c " OO " restart_ranking.results/reactions_discovery.log`

# Demonstrate LatticeDeformation
rm -rf complete_lattice_deformation.results

AMS_JOBNAME=complete_lattice_deformation $AMSBIN/reactions_discovery << eor
   MolecularDynamics
      Enabled Yes
      NumSimulations 2
      BuildSystem
         NumAtoms 50
         Density 0.4
         Molecule 
            SMILES C[CH2]  # a CH3CH2 radical
            MoleFraction 1
         End
         Molecule 
            SMILES [OH]   # an OH radical
            MoleFraction 3
         End
      End
      Type LatticeDeformation
      LatticeDeformation
         NumCycles 3
         Temperature 500
         MinVolumeFraction 0.2
         Period 200
      End
      BondOrders Method=Guess
   End

   NetworkExtraction Enabled=Yes UseCharges=No
   ProductRanking Enabled=Yes

   Engine ReaxFF
      ForceField CHON-2019.ff
      TaperBO Yes
   EndEngine
eor

#Test if ethanol is found:
echo ETHANOL: `grep -c " CCO " complete_lattice_deformation.results/reactions_discovery.log`
#Test if hydrogen peroxide is found:
echo HYDROGEN PEROXIDE: `grep -c " OO " complete_lattice_deformation.results/reactions_discovery.log`

# Read how many products there are in the output
nproducts=$($AMSBIN/amspython <<'eor'
from scm.base import KFFile
kf = KFFile("complete_nanoreactor.results/reactions_discovery.rkf")
products = kf.read("Categories", "Products")
print(len(products))
eor
)

# Print the paths created for the nproducts products
export AMS_JOBNAME=reaction_paths
$AMSBIN/reactions_discovery << eor
    SearchNetwork
        Enabled Yes
        ReactionNetwork complete_nanoreactor.results
$(for ((i=1; i<=nproducts; i++)); do
    echo "        ProductIndex $i"
done)
    End
eor

$AMSBIN/amspython << eor
from scm.base import KFFile

kf = KFFile("reaction_paths.results/reactions_discovery.rkf")
molindices = kf.read("ReactionNetworkResults", "MoleculeIndices")
npaths = len(molindices)
paths = [kf.read("ReactionNetworkResults", f"Path({i + 1})") for i in range(npaths)]

print("Cost info on ethanol and H2O2")
print("=============================")
for i, path in enumerate(paths):
    lines = path.split("\n")

    if lines[2].split()[-1] == "CCO":
        cost = float(lines[1].split()[-2])
        print(f"Cost Ethanol is {cost:10.5f} kcal/mol")

    if lines[2].split()[-1] == "OO":
        cost = float(lines[1].split()[-2])
        print(f"Cost H2O2 is {cost:10.5f} kcal/mol")
eor