{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "883f9446-cda0-4cf4-9b33-0b15ba825db9",
   "metadata": {},
   "source": [
    "# Example of Use\n",
    "\n",
    "These commands were used to generate the output data in [LAMMPS Simulation Data of Alchemical Processes](https://data.nist.gov/od/id/mds2-3637)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46d0cd27",
   "metadata": {},
   "source": [
    "## Import Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "451f49c8-8e17-4adf-9d09-21235461701f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "import generate_alchemical_lammps_inputs as gali\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bacb2eb",
   "metadata": {},
   "source": [
    "## Lennard-Jones Dimer\n",
    "\n",
    "A dimer made of two Lennard-Jones beads is solvated at T*=0.75. Here we generate a section of the LAMMPS input file to compute the solvation free energy of the dimer using MBAR. Note that this output will also contain the information to calculate the solvation free energy with thermodynamic integration also."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "78c98d7f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "# Variables and System Condition\n",
      "variable TK equal 0.7\n",
      "variable freq equal 100\n",
      "variable runtime_equil equal 100000\n",
      "variable runtime_prod equal 500000\n",
      "variable pinst equal pres\n",
      "variable tinst equal tem\n",
      "variable vinst equal vo\n",
      "variable pe equal p\n",
      "variable deltacdm equal 0.01 # delta used in central different method for derivative in T\n",
      "variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T\n",
      "variable nblocks equal \n",
      "variable lambdas vector [1.0,0.8,0.6,0.4,0.2,0.0\n",
      "\n",
      "thermo ${freq\n",
      "\n",
      "# Set-up Loo\n",
      "variable runid loop 1 ${nblocks} pa\n",
      "    label runloop\n",
      "\n",
      "    # Adjust param for the box and equilibrat\n",
      "    variable param equal v_lambdas[v_runid\n",
      "    if \"${runid} == 1\" then \n",
      "        \"jump SELF skipequil\n",
      "    variable ind equal v_runid-\n",
      "    variable param0 equal v_lambdas[v_ind\n",
      "    variable paramramp equal ramp(v_param0,v_param\n",
      "    fix ADAPT all adapt/fep 1 \n",
      "        pair lj/cut/soft lambda 1 2  v_paramramp\n",
      "    thermo_style custom step v_paramramp  temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Ram\n",
      "    unfix ADAP\n",
      "    fix ADAPT all adapt/fep ${freq} \n",
      "        pair lj/cut/soft lambda 1 2  v_param\n",
      "    thermo_style custom step v_param  temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Equi\n",
      "\n",
      "    label skipequil\n",
      "    write_data files/lj-cut-soft-lambda_${param}.dat\n",
      "\n",
      "    # Initialize compute\n",
      "    thermo_style custom step v_param temp press pe evdwl enthalp\n",
      "    compute FEPdb all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1 2  v_deltacdm2 \n",
      "        volume ye\n",
      "    compute FEPdf all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1 2  v_deltacdm \n",
      "        volume ye\n",
      "    fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_tinst v_pinst v_vinst v_pe \n",
      "        c_FEPdb[1] c_FEPdf[1] file files/ti_lj-cut-soft-lambda_${param}.txt\n",
      "    variable delta0 equal v_lambdas[1]-v_para\n",
      "    compute FEP000 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1 2  v_delta0 \n",
      "        volume ye\n",
      "    variable param000 equal v_param+v_delta\n",
      "    fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_pe \n",
      "        c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_lj-cut-soft-lambda_${param}_${param000}.txt\n",
      "    variable delta1 equal v_lambdas[2]-v_para\n",
      "    compute FEP001 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1 2  v_delta1 \n",
      "        volume ye\n",
      "    variable param001 equal v_param+v_delta\n",
      "    fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_pe \n",
      "        c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_lj-cut-soft-lambda_${param}_${param001}.txt\n",
      "    variable delta2 equal v_lambdas[3]-v_para\n",
      "    compute FEP002 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1 2  v_delta2 \n",
      "        volume ye\n",
      "    variable param002 equal v_param+v_delta\n",
      "    fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_pe \n",
      "        c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_lj-cut-soft-lambda_${param}_${param002}.txt\n",
      "    variable delta3 equal v_lambdas[4]-v_para\n",
      "    compute FEP003 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1 2  v_delta3 \n",
      "        volume ye\n",
      "    variable param003 equal v_param+v_delta\n",
      "    fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_pe \n",
      "        c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_lj-cut-soft-lambda_${param}_${param003}.txt\n",
      "    variable delta4 equal v_lambdas[5]-v_para\n",
      "    compute FEP004 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1 2  v_delta4 \n",
      "        volume ye\n",
      "    variable param004 equal v_param+v_delta\n",
      "    fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_pe \n",
      "        c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_lj-cut-soft-lambda_${param}_${param004}.txt\n",
      "    variable delta5 equal v_lambdas[6]-v_para\n",
      "    compute FEP005 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1 2  v_delta5 \n",
      "        volume ye\n",
      "    variable param005 equal v_param+v_delta\n",
      "    fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_pe \n",
      "        c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_lj-cut-soft-lambda_${param}_${param005}.txt\n",
      "  #  dump TRAJ all custom ${freq} files/dump_lj-cut-soft-lambda_${param}.lammpstrj id mol type element xu yu z\n",
      "\n",
      "    run ${runtime_prod}\n",
      "    uncompute FEPd\n",
      "    uncompute FEPd\n",
      "    if \"${runid} != 1\" then \n",
      "        \"unfix ADAPT\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    unfix FEPou\n",
      "  #  undump TRA\n",
      "\n",
      "    next runi\n",
      "    jump SELF runloop\n",
      "thermo_style custom v_param temp press pe evdwl enthalp\n",
      "write_data final.data nocoef\n"
     ]
    }
   ],
   "source": [
    "# Simple input where lambda pair changes\n",
    "output = gali.generate_mbar_input(\n",
    "    [1.0, 0.0], # The parameter will change from 1 to 0\n",
    "    -0.2, # The parameter will be reduced by 0.2 between each step, so there will be 5 windows\n",
    "    [\n",
    "        [ # This list represents an adaptation, more than one may be listed (e.g., change charge and nonbonded at the same time although that isn't recommended)\n",
    "            \"pair\", # May be \"pair\" or \"atom\", where charges can only be changed as per atom properties\n",
    "            (\n",
    "                \"lj/cut/soft\", # Pair style as defined in lammps\n",
    "                \"lambda\", # Define the changing parameter to be \"lambda\", an option for this pair style in the LAMMPS documentation (https://docs.lammps.org/fix_adapt_fep.html)\n",
    "                1, # Solvent bead type (or range)\n",
    "                2  # Sulute bead type (or range)\n",
    "            )\n",
    "        ]\n",
    "    ], \n",
    "    0.75, # Temperature\n",
    "#    output_file=\"your/output/here\", # Optional file to write to\n",
    "    n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)\n",
    "    n_run_prod_steps=5000000,  # Duration of production run for each step (i.e., window)\n",
    ")\n",
    "for line in output:\n",
    "    print(line[:-2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bdcf449",
   "metadata": {},
   "source": [
    "## Solvated Benzene"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fd0ed42",
   "metadata": {},
   "source": [
    "### Step 1: Scale Charges of Benzene to Zero\n",
    "\n",
    "The package LAMMPS cannot remove the Coulombic interactions between some atoms (e.g., benzene and water) while allowing other to remain (e.g., H in benzene with C in benzene). We must then scale the charges of the entire molecule to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "967ec82b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "# Variables and System Condition\n",
      "variable TK equal 30\n",
      "variable freq equal 100\n",
      "variable runtime_equil equal 100000\n",
      "variable runtime_prod equal 500000\n",
      "variable pinst equal pres\n",
      "variable tinst equal tem\n",
      "variable vinst equal vo\n",
      "variable pe equal p\n",
      "variable deltacdm equal 0.01 # delta used in central different method for derivative in T\n",
      "variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T\n",
      "variable nblocks equal \n",
      "variable lambdas vector [1.0,0.8,0.6,0.4,0.2,0.0\n",
      "variable deltacdmA0 equal -0.115*v_deltacd\n",
      "variable deltacdm2A0 equal 0.115*v_deltacd\n",
      "variable deltacdmA1 equal 0.115*v_deltacd\n",
      "variable deltacdm2A1 equal -0.115*v_deltacd\n",
      "\n",
      "thermo ${freq\n",
      "\n",
      "# Set-up Loo\n",
      "variable runid loop 1 ${nblocks} pa\n",
      "    label runloop\n",
      "\n",
      "    # Adjust param for the box and equilibrat\n",
      "    variable param equal v_lambdas[v_runid\n",
      "    if \"${runid} == 1\" then \n",
      "        \"jump SELF skipequil\n",
      "    variable ind equal v_runid-\n",
      "    variable param0 equal v_lambdas[v_ind\n",
      "    variable paramramp equal ramp(v_param0,v_param\n",
      "    variable paramrampA0 equal -0.115*ramp(v_param0,v_param\n",
      "    variable paramA0 equal -0.115*v_para\n",
      "    variable paramrampA1 equal 0.115*ramp(v_param0,v_param\n",
      "    variable paramA1 equal 0.115*v_para\n",
      "    fix ADAPT all adapt/fep 1 \n",
      "        atom charge 3  v_paramrampA0 \n",
      "        atom charge 4  v_paramrampA1\n",
      "    thermo_style custom step v_paramramp v_paramrampA0 v_paramrampA1 temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Ram\n",
      "    unfix ADAP\n",
      "    fix ADAPT all adapt/fep ${freq} \n",
      "        atom charge 3  v_paramA0 \n",
      "        atom charge 4  v_paramA1\n",
      "    thermo_style custom step v_param v_paramA0 v_paramA1 temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Equi\n",
      "\n",
      "    label skipequil\n",
      "    write_data files/charge_${param}.dat\n",
      "\n",
      "    # Initialize compute\n",
      "    thermo_style custom step v_param temp press pe evdwl enthalp\n",
      "    compute FEPdb all fep ${TK} \n",
      "        atom charge 3  v_deltacdm2A0 \n",
      "        atom charge 4  v_deltacdm2A1 \n",
      "        volume ye\n",
      "    compute FEPdf all fep ${TK} \n",
      "        atom charge 3  v_deltacdmA0 \n",
      "        atom charge 4  v_deltacdmA1 \n",
      "        volume ye\n",
      "    fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_tinst v_pinst v_vinst v_pe \n",
      "        c_FEPdb[1] c_FEPdf[1] file files/ti_charge_${param}.txt\n",
      "    variable delta0 equal v_lambdas[1]-v_para\n",
      "    variable deltaA00 equal -0.115*v_delta\n",
      "    variable deltaA01 equal 0.115*v_delta\n",
      "    compute FEP000 all fep ${TK} \n",
      "        atom charge 3  v_deltaA00 \n",
      "        atom charge 4  v_deltaA01 \n",
      "        volume ye\n",
      "    variable param000 equal v_param+v_delta\n",
      "    fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_pe \n",
      "        c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_charge_${param}_${param000}.txt\n",
      "    variable delta1 equal v_lambdas[2]-v_para\n",
      "    variable deltaA10 equal -0.115*v_delta\n",
      "    variable deltaA11 equal 0.115*v_delta\n",
      "    compute FEP001 all fep ${TK} \n",
      "        atom charge 3  v_deltaA10 \n",
      "        atom charge 4  v_deltaA11 \n",
      "        volume ye\n",
      "    variable param001 equal v_param+v_delta\n",
      "    fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_pe \n",
      "        c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_charge_${param}_${param001}.txt\n",
      "    variable delta2 equal v_lambdas[3]-v_para\n",
      "    variable deltaA20 equal -0.115*v_delta\n",
      "    variable deltaA21 equal 0.115*v_delta\n",
      "    compute FEP002 all fep ${TK} \n",
      "        atom charge 3  v_deltaA20 \n",
      "        atom charge 4  v_deltaA21 \n",
      "        volume ye\n",
      "    variable param002 equal v_param+v_delta\n",
      "    fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_pe \n",
      "        c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_charge_${param}_${param002}.txt\n",
      "    variable delta3 equal v_lambdas[4]-v_para\n",
      "    variable deltaA30 equal -0.115*v_delta\n",
      "    variable deltaA31 equal 0.115*v_delta\n",
      "    compute FEP003 all fep ${TK} \n",
      "        atom charge 3  v_deltaA30 \n",
      "        atom charge 4  v_deltaA31 \n",
      "        volume ye\n",
      "    variable param003 equal v_param+v_delta\n",
      "    fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_pe \n",
      "        c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_charge_${param}_${param003}.txt\n",
      "    variable delta4 equal v_lambdas[5]-v_para\n",
      "    variable deltaA40 equal -0.115*v_delta\n",
      "    variable deltaA41 equal 0.115*v_delta\n",
      "    compute FEP004 all fep ${TK} \n",
      "        atom charge 3  v_deltaA40 \n",
      "        atom charge 4  v_deltaA41 \n",
      "        volume ye\n",
      "    variable param004 equal v_param+v_delta\n",
      "    fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_pe \n",
      "        c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_charge_${param}_${param004}.txt\n",
      "    variable delta5 equal v_lambdas[6]-v_para\n",
      "    variable deltaA50 equal -0.115*v_delta\n",
      "    variable deltaA51 equal 0.115*v_delta\n",
      "    compute FEP005 all fep ${TK} \n",
      "        atom charge 3  v_deltaA50 \n",
      "        atom charge 4  v_deltaA51 \n",
      "        volume ye\n",
      "    variable param005 equal v_param+v_delta\n",
      "    fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_pe \n",
      "        c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_charge_${param}_${param005}.txt\n",
      "  #  dump TRAJ all custom ${freq} files/dump_charge_${param}.lammpstrj id mol type element xu yu z\n",
      "\n",
      "    run ${runtime_prod}\n",
      "    uncompute FEPd\n",
      "    uncompute FEPd\n",
      "    if \"${runid} != 1\" then \n",
      "        \"unfix ADAPT\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    unfix FEPou\n",
      "  #  undump TRA\n",
      "\n",
      "    next runi\n",
      "    jump SELF runloop\n",
      "thermo_style custom v_param temp press pe evdwl enthalp\n",
      "write_data final.data nocoef\n"
     ]
    }
   ],
   "source": [
    "# Input where charges change\n",
    "output = gali.generate_mbar_input(\n",
    "    [1.0, 0.0], # Scale charges from full strength to none\n",
    "    -0.2, # Decrease in steps of 0.2 (5 windows)\n",
    "    [\n",
    "        [\"atom\", (\"charge\", 3, -0.115000)], # Change the atom property, charge, of atom type 3, with a full strength value of -0.115000\n",
    "        [\"atom\", (\"charge\", 4, 0.115000)], # Change another atom property...\n",
    "    ],\n",
    "    300, # Temperature in K\n",
    "#    output_file=\"your/output/here\", # Optional file to write to\n",
    "    n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)\n",
    "    n_run_prod_steps=5000000,  # Duration of production run for each step (i.e., window)\n",
    ")\n",
    "for line in output:\n",
    "    print(line[:-2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baa6cfd5",
   "metadata": {},
   "source": [
    "### Step 2a: Remove Non-bonded Interaction between Benzene and Water without Trajectory\n",
    "\n",
    "Because saving the trajectories for a free energy calculation can be quite large, this is the preferred method. However, if later you determine that you need an additional value of lambda, you'll need to rerun the entire calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "780ed25b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "# Variables and System Condition\n",
      "variable TK equal 30\n",
      "variable freq equal 100\n",
      "variable runtime_equil equal 100000\n",
      "variable runtime_prod equal 500000\n",
      "variable pinst equal pres\n",
      "variable tinst equal tem\n",
      "variable vinst equal vo\n",
      "variable pe equal p\n",
      "variable deltacdm equal 0.01 # delta used in central different method for derivative in T\n",
      "variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T\n",
      "variable nblocks equal 1\n",
      "variable lambdas vector [1.0,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0\n",
      "\n",
      "# Set Previous Chang\n",
      "variable param2 equal \n",
      "variable param2A0 equal -0.115*v_param\n",
      "variable param2A1 equal 0.115*v_param\n",
      "fix ADAPT2 all adapt/fep 1 \n",
      "    atom charge 3  v_param2A0 \n",
      "    atom charge 4  v_param2A1\n",
      "\n",
      "thermo ${freq\n",
      "\n",
      "# Set-up Loo\n",
      "variable runid loop 1 ${nblocks} pa\n",
      "    label runloop\n",
      "\n",
      "    # Adjust param for the box and equilibrat\n",
      "    variable param equal v_lambdas[v_runid\n",
      "    if \"${runid} == 1\" then \n",
      "        \"jump SELF skipequil\n",
      "    variable ind equal v_runid-\n",
      "    variable param0 equal v_lambdas[v_ind\n",
      "    variable paramramp equal ramp(v_param0,v_param\n",
      "    fix ADAPT all adapt/fep 1 \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_paramramp\n",
      "    thermo_style custom step v_paramramp  temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Ram\n",
      "    unfix ADAP\n",
      "    fix ADAPT all adapt/fep ${freq} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_param\n",
      "    thermo_style custom step v_param  temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Equi\n",
      "\n",
      "    label skipequil\n",
      "    write_data files/lj-cut-soft-lambda_${param}.dat\n",
      "\n",
      "    # Initialize compute\n",
      "    thermo_style custom step v_param temp press pe evdwl enthalp\n",
      "    compute FEPdb all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_deltacdm2 \n",
      "        volume ye\n",
      "    compute FEPdf all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_deltacdm \n",
      "        volume ye\n",
      "    fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_param2 v_tinst v_pinst v_vinst v_pe \n",
      "        c_FEPdb[1] c_FEPdf[1] file files/ti_lj-cut-soft-lambda_${param}_charge--0.0_0.txt\n",
      "    variable delta0 equal v_lambdas[1]-v_para\n",
      "    compute FEP000 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta0 \n",
      "        volume ye\n",
      "    variable param000 equal v_param+v_delta\n",
      "    fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_param2 v_pe \n",
      "        c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_lj-cut-soft-lambda_${param}_${param000}_charge--0.0.txt\n",
      "    variable delta1 equal v_lambdas[2]-v_para\n",
      "    compute FEP001 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta1 \n",
      "        volume ye\n",
      "    variable param001 equal v_param+v_delta\n",
      "    fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_param2 v_pe \n",
      "        c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_lj-cut-soft-lambda_${param}_${param001}_charge--0.0.txt\n",
      "    variable delta2 equal v_lambdas[3]-v_para\n",
      "    compute FEP002 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta2 \n",
      "        volume ye\n",
      "    variable param002 equal v_param+v_delta\n",
      "    fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_param2 v_pe \n",
      "        c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_lj-cut-soft-lambda_${param}_${param002}_charge--0.0.txt\n",
      "    variable delta3 equal v_lambdas[4]-v_para\n",
      "    compute FEP003 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta3 \n",
      "        volume ye\n",
      "    variable param003 equal v_param+v_delta\n",
      "    fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_param2 v_pe \n",
      "        c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_lj-cut-soft-lambda_${param}_${param003}_charge--0.0.txt\n",
      "    variable delta4 equal v_lambdas[5]-v_para\n",
      "    compute FEP004 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta4 \n",
      "        volume ye\n",
      "    variable param004 equal v_param+v_delta\n",
      "    fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_param2 v_pe \n",
      "        c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_lj-cut-soft-lambda_${param}_${param004}_charge--0.0.txt\n",
      "    variable delta5 equal v_lambdas[6]-v_para\n",
      "    compute FEP005 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta5 \n",
      "        volume ye\n",
      "    variable param005 equal v_param+v_delta\n",
      "    fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_param2 v_pe \n",
      "        c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_lj-cut-soft-lambda_${param}_${param005}_charge--0.0.txt\n",
      "    variable delta6 equal v_lambdas[7]-v_para\n",
      "    compute FEP006 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta6 \n",
      "        volume ye\n",
      "    variable param006 equal v_param+v_delta\n",
      "    fix FEPout006 all ave/time ${freq} 1 ${freq} v_param v_param006 v_param2 v_pe \n",
      "        c_FEP006[1] c_FEP006[2] c_FEP006[3] file files/mbar_lj-cut-soft-lambda_${param}_${param006}_charge--0.0.txt\n",
      "    variable delta7 equal v_lambdas[8]-v_para\n",
      "    compute FEP007 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta7 \n",
      "        volume ye\n",
      "    variable param007 equal v_param+v_delta\n",
      "    fix FEPout007 all ave/time ${freq} 1 ${freq} v_param v_param007 v_param2 v_pe \n",
      "        c_FEP007[1] c_FEP007[2] c_FEP007[3] file files/mbar_lj-cut-soft-lambda_${param}_${param007}_charge--0.0.txt\n",
      "    variable delta8 equal v_lambdas[9]-v_para\n",
      "    compute FEP008 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta8 \n",
      "        volume ye\n",
      "    variable param008 equal v_param+v_delta\n",
      "    fix FEPout008 all ave/time ${freq} 1 ${freq} v_param v_param008 v_param2 v_pe \n",
      "        c_FEP008[1] c_FEP008[2] c_FEP008[3] file files/mbar_lj-cut-soft-lambda_${param}_${param008}_charge--0.0.txt\n",
      "    variable delta9 equal v_lambdas[10]-v_para\n",
      "    compute FEP009 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta9 \n",
      "        volume ye\n",
      "    variable param009 equal v_param+v_delta\n",
      "    fix FEPout009 all ave/time ${freq} 1 ${freq} v_param v_param009 v_param2 v_pe \n",
      "        c_FEP009[1] c_FEP009[2] c_FEP009[3] file files/mbar_lj-cut-soft-lambda_${param}_${param009}_charge--0.0.txt\n",
      "    variable delta10 equal v_lambdas[11]-v_para\n",
      "    compute FEP010 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta10 \n",
      "        volume ye\n",
      "    variable param010 equal v_param+v_delta1\n",
      "    fix FEPout010 all ave/time ${freq} 1 ${freq} v_param v_param010 v_param2 v_pe \n",
      "        c_FEP010[1] c_FEP010[2] c_FEP010[3] file files/mbar_lj-cut-soft-lambda_${param}_${param010}_charge--0.0.txt\n",
      "    variable delta11 equal v_lambdas[12]-v_para\n",
      "    compute FEP011 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta11 \n",
      "        volume ye\n",
      "    variable param011 equal v_param+v_delta1\n",
      "    fix FEPout011 all ave/time ${freq} 1 ${freq} v_param v_param011 v_param2 v_pe \n",
      "        c_FEP011[1] c_FEP011[2] c_FEP011[3] file files/mbar_lj-cut-soft-lambda_${param}_${param011}_charge--0.0.txt\n",
      "    variable delta12 equal v_lambdas[13]-v_para\n",
      "    compute FEP012 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta12 \n",
      "        volume ye\n",
      "    variable param012 equal v_param+v_delta1\n",
      "    fix FEPout012 all ave/time ${freq} 1 ${freq} v_param v_param012 v_param2 v_pe \n",
      "        c_FEP012[1] c_FEP012[2] c_FEP012[3] file files/mbar_lj-cut-soft-lambda_${param}_${param012}_charge--0.0.txt\n",
      "    variable delta13 equal v_lambdas[14]-v_para\n",
      "    compute FEP013 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta13 \n",
      "        volume ye\n",
      "    variable param013 equal v_param+v_delta1\n",
      "    fix FEPout013 all ave/time ${freq} 1 ${freq} v_param v_param013 v_param2 v_pe \n",
      "        c_FEP013[1] c_FEP013[2] c_FEP013[3] file files/mbar_lj-cut-soft-lambda_${param}_${param013}_charge--0.0.txt\n",
      "    variable delta14 equal v_lambdas[15]-v_para\n",
      "    compute FEP014 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta14 \n",
      "        volume ye\n",
      "    variable param014 equal v_param+v_delta1\n",
      "    fix FEPout014 all ave/time ${freq} 1 ${freq} v_param v_param014 v_param2 v_pe \n",
      "        c_FEP014[1] c_FEP014[2] c_FEP014[3] file files/mbar_lj-cut-soft-lambda_${param}_${param014}_charge--0.0.txt\n",
      "    variable delta15 equal v_lambdas[16]-v_para\n",
      "    compute FEP015 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta15 \n",
      "        volume ye\n",
      "    variable param015 equal v_param+v_delta1\n",
      "    fix FEPout015 all ave/time ${freq} 1 ${freq} v_param v_param015 v_param2 v_pe \n",
      "        c_FEP015[1] c_FEP015[2] c_FEP015[3] file files/mbar_lj-cut-soft-lambda_${param}_${param015}_charge--0.0.txt\n",
      "  #  dump TRAJ all custom ${freq} files/dump_lj-cut-soft-lambda_${param}_charge--0.0_0.lammpstrj id mol type element xu yu z\n",
      "\n",
      "    run ${runtime_prod}\n",
      "    uncompute FEPd\n",
      "    uncompute FEPd\n",
      "    if \"${runid} != 1\" then \n",
      "        \"unfix ADAPT\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    unfix FEPou\n",
      "  #  undump TRA\n",
      "\n",
      "    next runi\n",
      "    jump SELF runloop\n",
      "thermo_style custom v_param temp press pe evdwl enthalp\n",
      "write_data final.data nocoef\n"
     ]
    }
   ],
   "source": [
    "# Input where pair changes after charges change\n",
    "output = gali.generate_mbar_input(\n",
    "    [1.0, 0.0], \n",
    "    -0.0625, # Remove within 16 steps\n",
    "    [\n",
    "        [\n",
    "            \"pair\", # change a pair style (in the previous step we changed an atom property)\n",
    "            (\n",
    "                \"lj/cut/soft\",\n",
    "                \"lambda\", \n",
    "                \"1*2\", # Range of atom types for water (in format required by LAMMPS)\n",
    "                \"3*4\" # Range of atom types for the solute, benzene (in format required by LAMMPS)\n",
    "            )\n",
    "        ]\n",
    "    ],\n",
    "    300,\n",
    "    parameter_array=np.array([1,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0]), # Overwrite parameter_change with custom lambda spacing\n",
    "    fix_adapt_changes2 = [ # Keep the change in charges done in the previous step, so defined again here\n",
    "        [\"atom\", (\"charge\", 3, -0.115000)],\n",
    "        [\"atom\", (\"charge\", 4, 0.115000)],   \n",
    "    ],\n",
    "    parameter2_value=0, # Define the final scaling factor for the parameter change in step 1\n",
    "#    output_file=\"your/output/here\", # Optional file to write to\n",
    "    n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)\n",
    "    n_run_prod_steps=5000000,  # Duration of production run for each step (i.e., window)\n",
    ")\n",
    "for line in output:\n",
    "    print(line[:-2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de30ae76",
   "metadata": {},
   "source": [
    "### Step 2b: Remove Non-bonded Interaction between Benzene and Water with Trajectory\n",
    "\n",
    "Generate the input file to write the trajectories needed to calculate the free energy of solvation. In the next step, we generate the input to calculate the values needed for MBAR FOR ONE WINDOW using rerun and the trajectories. This might be desirable if you don't know the number of lambda windows or spacing that your system requires."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "66e7ecb7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "# Variables and System Condition\n",
      "variable TK equal 30\n",
      "variable freq equal 100\n",
      "variable runtime_equil equal 100000\n",
      "variable runtime_prod equal 500000\n",
      "variable pinst equal pres\n",
      "variable tinst equal tem\n",
      "variable vinst equal vo\n",
      "variable pe equal p\n",
      "variable deltacdm equal 0.01 # delta used in central different method for derivative in T\n",
      "variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T\n",
      "variable nblocks equal 1\n",
      "variable lambdas vector [1.0,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0\n",
      "\n",
      "# Set Previous Chang\n",
      "variable param2 equal \n",
      "variable param2A0 equal -0.115*v_param\n",
      "variable param2A1 equal 0.115*v_param\n",
      "fix ADAPT2 all adapt/fep 1 \n",
      "    atom charge 3  v_param2A0 \n",
      "    atom charge 4  v_param2A1\n",
      "\n",
      "thermo ${freq\n",
      "\n",
      "# Set-up Loo\n",
      "variable runid loop 1 ${nblocks} pa\n",
      "    label runloop\n",
      "\n",
      "    # Adjust param for the box and equilibrat\n",
      "    variable param equal v_lambdas[v_runid\n",
      "    if \"${runid} == 1\" then \n",
      "        \"jump SELF skipequil\n",
      "    variable ind equal v_runid-\n",
      "    variable param0 equal v_lambdas[v_ind\n",
      "    variable paramramp equal ramp(v_param0,v_param\n",
      "    fix ADAPT all adapt/fep 1 \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_paramramp\n",
      "    thermo_style custom step v_paramramp  temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Ram\n",
      "    unfix ADAP\n",
      "    fix ADAPT all adapt/fep ${freq} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_param\n",
      "    thermo_style custom step v_param  temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Equi\n",
      "\n",
      "    label skipequil\n",
      "    write_data files/lj-cut-soft-lambda_${param}.dat\n",
      "\n",
      "    # Initialize compute\n",
      "    thermo_style custom step v_param temp press pe evdwl enthalp\n",
      "    compute FEPdb all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_deltacdm2 \n",
      "        volume ye\n",
      "    compute FEPdf all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_deltacdm \n",
      "        volume ye\n",
      "    fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_param2 v_tinst v_pinst v_vinst v_pe \n",
      "        c_FEPdb[1] c_FEPdf[1] file files/ti_lj-cut-soft-lambda_${param}_charge--0.0_0.txt\n",
      "    variable delta0 equal v_lambdas[1]-v_para\n",
      "    compute FEP000 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta0 \n",
      "        volume ye\n",
      "    variable param000 equal v_param+v_delta\n",
      "    fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_param2 v_pe \n",
      "        c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_lj-cut-soft-lambda_${param}_${param000}_charge--0.0.txt\n",
      "    variable delta1 equal v_lambdas[2]-v_para\n",
      "    compute FEP001 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta1 \n",
      "        volume ye\n",
      "    variable param001 equal v_param+v_delta\n",
      "    fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_param2 v_pe \n",
      "        c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_lj-cut-soft-lambda_${param}_${param001}_charge--0.0.txt\n",
      "    variable delta2 equal v_lambdas[3]-v_para\n",
      "    compute FEP002 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta2 \n",
      "        volume ye\n",
      "    variable param002 equal v_param+v_delta\n",
      "    fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_param2 v_pe \n",
      "        c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_lj-cut-soft-lambda_${param}_${param002}_charge--0.0.txt\n",
      "    variable delta3 equal v_lambdas[4]-v_para\n",
      "    compute FEP003 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta3 \n",
      "        volume ye\n",
      "    variable param003 equal v_param+v_delta\n",
      "    fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_param2 v_pe \n",
      "        c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_lj-cut-soft-lambda_${param}_${param003}_charge--0.0.txt\n",
      "    variable delta4 equal v_lambdas[5]-v_para\n",
      "    compute FEP004 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta4 \n",
      "        volume ye\n",
      "    variable param004 equal v_param+v_delta\n",
      "    fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_param2 v_pe \n",
      "        c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_lj-cut-soft-lambda_${param}_${param004}_charge--0.0.txt\n",
      "    variable delta5 equal v_lambdas[6]-v_para\n",
      "    compute FEP005 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta5 \n",
      "        volume ye\n",
      "    variable param005 equal v_param+v_delta\n",
      "    fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_param2 v_pe \n",
      "        c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_lj-cut-soft-lambda_${param}_${param005}_charge--0.0.txt\n",
      "    variable delta6 equal v_lambdas[7]-v_para\n",
      "    compute FEP006 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta6 \n",
      "        volume ye\n",
      "    variable param006 equal v_param+v_delta\n",
      "    fix FEPout006 all ave/time ${freq} 1 ${freq} v_param v_param006 v_param2 v_pe \n",
      "        c_FEP006[1] c_FEP006[2] c_FEP006[3] file files/mbar_lj-cut-soft-lambda_${param}_${param006}_charge--0.0.txt\n",
      "    variable delta7 equal v_lambdas[8]-v_para\n",
      "    compute FEP007 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta7 \n",
      "        volume ye\n",
      "    variable param007 equal v_param+v_delta\n",
      "    fix FEPout007 all ave/time ${freq} 1 ${freq} v_param v_param007 v_param2 v_pe \n",
      "        c_FEP007[1] c_FEP007[2] c_FEP007[3] file files/mbar_lj-cut-soft-lambda_${param}_${param007}_charge--0.0.txt\n",
      "    variable delta8 equal v_lambdas[9]-v_para\n",
      "    compute FEP008 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta8 \n",
      "        volume ye\n",
      "    variable param008 equal v_param+v_delta\n",
      "    fix FEPout008 all ave/time ${freq} 1 ${freq} v_param v_param008 v_param2 v_pe \n",
      "        c_FEP008[1] c_FEP008[2] c_FEP008[3] file files/mbar_lj-cut-soft-lambda_${param}_${param008}_charge--0.0.txt\n",
      "    variable delta9 equal v_lambdas[10]-v_para\n",
      "    compute FEP009 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta9 \n",
      "        volume ye\n",
      "    variable param009 equal v_param+v_delta\n",
      "    fix FEPout009 all ave/time ${freq} 1 ${freq} v_param v_param009 v_param2 v_pe \n",
      "        c_FEP009[1] c_FEP009[2] c_FEP009[3] file files/mbar_lj-cut-soft-lambda_${param}_${param009}_charge--0.0.txt\n",
      "    variable delta10 equal v_lambdas[11]-v_para\n",
      "    compute FEP010 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta10 \n",
      "        volume ye\n",
      "    variable param010 equal v_param+v_delta1\n",
      "    fix FEPout010 all ave/time ${freq} 1 ${freq} v_param v_param010 v_param2 v_pe \n",
      "        c_FEP010[1] c_FEP010[2] c_FEP010[3] file files/mbar_lj-cut-soft-lambda_${param}_${param010}_charge--0.0.txt\n",
      "    variable delta11 equal v_lambdas[12]-v_para\n",
      "    compute FEP011 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta11 \n",
      "        volume ye\n",
      "    variable param011 equal v_param+v_delta1\n",
      "    fix FEPout011 all ave/time ${freq} 1 ${freq} v_param v_param011 v_param2 v_pe \n",
      "        c_FEP011[1] c_FEP011[2] c_FEP011[3] file files/mbar_lj-cut-soft-lambda_${param}_${param011}_charge--0.0.txt\n",
      "    variable delta12 equal v_lambdas[13]-v_para\n",
      "    compute FEP012 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta12 \n",
      "        volume ye\n",
      "    variable param012 equal v_param+v_delta1\n",
      "    fix FEPout012 all ave/time ${freq} 1 ${freq} v_param v_param012 v_param2 v_pe \n",
      "        c_FEP012[1] c_FEP012[2] c_FEP012[3] file files/mbar_lj-cut-soft-lambda_${param}_${param012}_charge--0.0.txt\n",
      "    variable delta13 equal v_lambdas[14]-v_para\n",
      "    compute FEP013 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta13 \n",
      "        volume ye\n",
      "    variable param013 equal v_param+v_delta1\n",
      "    fix FEPout013 all ave/time ${freq} 1 ${freq} v_param v_param013 v_param2 v_pe \n",
      "        c_FEP013[1] c_FEP013[2] c_FEP013[3] file files/mbar_lj-cut-soft-lambda_${param}_${param013}_charge--0.0.txt\n",
      "    variable delta14 equal v_lambdas[15]-v_para\n",
      "    compute FEP014 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta14 \n",
      "        volume ye\n",
      "    variable param014 equal v_param+v_delta1\n",
      "    fix FEPout014 all ave/time ${freq} 1 ${freq} v_param v_param014 v_param2 v_pe \n",
      "        c_FEP014[1] c_FEP014[2] c_FEP014[3] file files/mbar_lj-cut-soft-lambda_${param}_${param014}_charge--0.0.txt\n",
      "    variable delta15 equal v_lambdas[16]-v_para\n",
      "    compute FEP015 all fep ${TK} \n",
      "        pair lj/cut/soft lambda 1*2 3*4  v_delta15 \n",
      "        volume ye\n",
      "    variable param015 equal v_param+v_delta1\n",
      "    fix FEPout015 all ave/time ${freq} 1 ${freq} v_param v_param015 v_param2 v_pe \n",
      "        c_FEP015[1] c_FEP015[2] c_FEP015[3] file files/mbar_lj-cut-soft-lambda_${param}_${param015}_charge--0.0.txt\n",
      "  #  dump TRAJ all custom ${freq} files/dump_lj-cut-soft-lambda_${param}_charge--0.0_0.lammpstrj id mol type element xu yu z\n",
      "\n",
      "    run ${runtime_prod}\n",
      "    uncompute FEPd\n",
      "    uncompute FEPd\n",
      "    if \"${runid} != 1\" then \n",
      "        \"unfix ADAPT\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    uncompute FEP01\n",
      "    unfix FEPout01\n",
      "    unfix FEPou\n",
      "  #  undump TRA\n",
      "\n",
      "    next runi\n",
      "    jump SELF runloop\n",
      "thermo_style custom v_param temp press pe evdwl enthalp\n",
      "write_data final.data nocoef\n"
     ]
    }
   ],
   "source": [
    "# Input where pair changes after charges change\n",
    "ouptut = gali.generate_traj_input(\n",
    "    [1.0, 0.0], \n",
    "    -0.0625, \n",
    "    [\n",
    "        [\"pair\", (\"lj/cut/soft\", \"lambda\", \"1*2\", \"3*4\")]\n",
    "    ],\n",
    "    300,\n",
    "    parameter_array=np.array([1,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0]), # Overwrite parameter_change with custom lambda spacing\n",
    "    fix_adapt_changes2 = [\n",
    "        [\"atom\", (\"charge\", 3, -0.115000)],\n",
    "        [\"atom\", (\"charge\", 4, 0.115000)],   \n",
    "    ],\n",
    "    parameter2_value=0,\n",
    "#    output_file=\"your/output/here\", # Optional file to write to\n",
    "    n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)\n",
    "    n_run_prod_steps=5000000,  # Duration of production run for each step (i.e., window)\n",
    ")\n",
    "for line in output:\n",
    "    print(line[:-2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4b6d7210",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "# Variables and System Condition\n",
      "variable TK equal 30\n",
      "variable freq equal 100\n",
      "variable pinst equal pres\n",
      "variable tinst equal tem\n",
      "variable vinst equal vo\n",
      "variable pe equal p\n",
      "variable param equal 0.\n",
      "variable lambdas vector [1.0,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0\n",
      "\n",
      "thermo ${freq\n",
      "\n",
      "# Set Previous Chang\n",
      "variable param2 equal \n",
      "variable param2A0 equal -0.115*v_param\n",
      "variable param2A1 equal 0.115*v_param\n",
      "fix ADAPT2 all adapt/fep 1 \n",
      "    atom charge 3  v_param2A0 \n",
      "    atom charge 4  v_param2A1\n",
      "\n",
      "# Initialize compute\n",
      "fix ADAPT all adapt/fep ${freq} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_param \n",
      "variable delta0 equal 0.\n",
      "compute FEP000 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta0 \n",
      "    volume ye\n",
      "variable param000 equal v_param+v_delta\n",
      "fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_param2 v_pe \n",
      "    c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_lj-cut-soft-lambda_${param}_${param000}_charge--0.0.txt\n",
      "variable delta1 equal 0.4\n",
      "compute FEP001 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta1 \n",
      "    volume ye\n",
      "variable param001 equal v_param+v_delta\n",
      "fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_param2 v_pe \n",
      "    c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_lj-cut-soft-lambda_${param}_${param001}_charge--0.0.txt\n",
      "variable delta2 equal 0.\n",
      "compute FEP002 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta2 \n",
      "    volume ye\n",
      "variable param002 equal v_param+v_delta\n",
      "fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_param2 v_pe \n",
      "    c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_lj-cut-soft-lambda_${param}_${param002}_charge--0.0.txt\n",
      "variable delta3 equal 0.3\n",
      "compute FEP003 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta3 \n",
      "    volume ye\n",
      "variable param003 equal v_param+v_delta\n",
      "fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_param2 v_pe \n",
      "    c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_lj-cut-soft-lambda_${param}_${param003}_charge--0.0.txt\n",
      "variable delta4 equal 0.\n",
      "compute FEP004 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta4 \n",
      "    volume ye\n",
      "variable param004 equal v_param+v_delta\n",
      "fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_param2 v_pe \n",
      "    c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_lj-cut-soft-lambda_${param}_${param004}_charge--0.0.txt\n",
      "variable delta5 equal 0.2\n",
      "compute FEP005 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta5 \n",
      "    volume ye\n",
      "variable param005 equal v_param+v_delta\n",
      "fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_param2 v_pe \n",
      "    c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_lj-cut-soft-lambda_${param}_${param005}_charge--0.0.txt\n",
      "variable delta6 equal 0.\n",
      "compute FEP006 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta6 \n",
      "    volume ye\n",
      "variable param006 equal v_param+v_delta\n",
      "fix FEPout006 all ave/time ${freq} 1 ${freq} v_param v_param006 v_param2 v_pe \n",
      "    c_FEP006[1] c_FEP006[2] c_FEP006[3] file files/mbar_lj-cut-soft-lambda_${param}_${param006}_charge--0.0.txt\n",
      "variable delta7 equal 0.1\n",
      "compute FEP007 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta7 \n",
      "    volume ye\n",
      "variable param007 equal v_param+v_delta\n",
      "fix FEPout007 all ave/time ${freq} 1 ${freq} v_param v_param007 v_param2 v_pe \n",
      "    c_FEP007[1] c_FEP007[2] c_FEP007[3] file files/mbar_lj-cut-soft-lambda_${param}_${param007}_charge--0.0.txt\n",
      "variable delta8 equal 0.\n",
      "compute FEP008 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta8 \n",
      "    volume ye\n",
      "variable param008 equal v_param+v_delta\n",
      "fix FEPout008 all ave/time ${freq} 1 ${freq} v_param v_param008 v_param2 v_pe \n",
      "    c_FEP008[1] c_FEP008[2] c_FEP008[3] file files/mbar_lj-cut-soft-lambda_${param}_${param008}_charge--0.0.txt\n",
      "variable delta9 equal 0.\n",
      "compute FEP009 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta9 \n",
      "    volume ye\n",
      "variable param009 equal v_param+v_delta\n",
      "fix FEPout009 all ave/time ${freq} 1 ${freq} v_param v_param009 v_param2 v_pe \n",
      "    c_FEP009[1] c_FEP009[2] c_FEP009[3] file files/mbar_lj-cut-soft-lambda_${param}_${param009}_charge--0.0.txt\n",
      "variable delta10 equal -0.\n",
      "compute FEP010 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta10 \n",
      "    volume ye\n",
      "variable param010 equal v_param+v_delta1\n",
      "fix FEPout010 all ave/time ${freq} 1 ${freq} v_param v_param010 v_param2 v_pe \n",
      "    c_FEP010[1] c_FEP010[2] c_FEP010[3] file files/mbar_lj-cut-soft-lambda_${param}_${param010}_charge--0.0.txt\n",
      "variable delta11 equal -0.\n",
      "compute FEP011 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta11 \n",
      "    volume ye\n",
      "variable param011 equal v_param+v_delta1\n",
      "fix FEPout011 all ave/time ${freq} 1 ${freq} v_param v_param011 v_param2 v_pe \n",
      "    c_FEP011[1] c_FEP011[2] c_FEP011[3] file files/mbar_lj-cut-soft-lambda_${param}_${param011}_charge--0.0.txt\n",
      "variable delta12 equal -0.\n",
      "compute FEP012 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta12 \n",
      "    volume ye\n",
      "variable param012 equal v_param+v_delta1\n",
      "fix FEPout012 all ave/time ${freq} 1 ${freq} v_param v_param012 v_param2 v_pe \n",
      "    c_FEP012[1] c_FEP012[2] c_FEP012[3] file files/mbar_lj-cut-soft-lambda_${param}_${param012}_charge--0.0.txt\n",
      "variable delta13 equal -0.\n",
      "compute FEP013 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta13 \n",
      "    volume ye\n",
      "variable param013 equal v_param+v_delta1\n",
      "fix FEPout013 all ave/time ${freq} 1 ${freq} v_param v_param013 v_param2 v_pe \n",
      "    c_FEP013[1] c_FEP013[2] c_FEP013[3] file files/mbar_lj-cut-soft-lambda_${param}_${param013}_charge--0.0.txt\n",
      "variable delta14 equal -0.4\n",
      "compute FEP014 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta14 \n",
      "    volume ye\n",
      "variable param014 equal v_param+v_delta1\n",
      "fix FEPout014 all ave/time ${freq} 1 ${freq} v_param v_param014 v_param2 v_pe \n",
      "    c_FEP014[1] c_FEP014[2] c_FEP014[3] file files/mbar_lj-cut-soft-lambda_${param}_${param014}_charge--0.0.txt\n",
      "variable delta15 equal -0.\n",
      "compute FEP015 all fep ${TK} \n",
      "    pair lj/cut/soft lambda 1*2 3*4  v_delta15 \n",
      "    volume ye\n",
      "variable param015 equal v_param+v_delta1\n",
      "fix FEPout015 all ave/time ${freq} 1 ${freq} v_param v_param015 v_param2 v_pe \n",
      "    c_FEP015[1] c_FEP015[2] c_FEP015[3] file files/mbar_lj-cut-soft-lambda_${param}_${param015}_charge--0.0.txt\n",
      "\n",
      "rerun files/dump_lj-cut-soft-lambda_${param}_charge--0.0_0.lammpstrj every ${freq} dump x y z\n"
     ]
    }
   ],
   "source": [
    "output = gali.generate_rerun_mbar(\n",
    "    0.5, # Value of lambda to make a rerun script for\n",
    "    [1.0, 0.0], \n",
    "    -0.0625, \n",
    "    [\n",
    "        [\"pair\", (\"lj/cut/soft\", \"lambda\", \"1*2\", \"3*4\")]\n",
    "    ],\n",
    "    300,\n",
    "    parameter_array=np.array([1,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0]),\n",
    "    fix_adapt_changes2 = [\n",
    "        [\"atom\", (\"charge\", 3, -0.115000)],\n",
    "        [\"atom\", (\"charge\", 4, 0.115000)],   \n",
    "    ],\n",
    "    parameter2_value=0,\n",
    "#    output_file=\"your/output/here\", # Optional file to write to\n",
    ")\n",
    "for line in output:\n",
    "    print(line[:-2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba8158d8",
   "metadata": {},
   "source": [
    "### Step 3: Restore the Charges of Benzene Atoms in Vacuum\n",
    "\n",
    "Because the previous two steps were missing the intramolecular charge interactions of benzene, and the end of the second system resulted in benzene being completely decoupled from water, we must now determine the free energy of returning the charges to complete the cycle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "562ce1e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "# Variables and System Condition\n",
      "variable TK equal 30\n",
      "variable freq equal 100\n",
      "variable runtime_equil equal 100000\n",
      "variable runtime_prod equal 500000\n",
      "variable pinst equal pres\n",
      "variable tinst equal tem\n",
      "variable vinst equal vo\n",
      "variable pe equal p\n",
      "variable deltacdm equal 0.01 # delta used in central different method for derivative in T\n",
      "variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T\n",
      "variable nblocks equal \n",
      "variable lambdas vector [0.0,0.2,0.4,0.6,0.8,1.0\n",
      "variable deltacdmA0 equal -0.115*v_deltacd\n",
      "variable deltacdm2A0 equal 0.115*v_deltacd\n",
      "variable deltacdmA1 equal 0.115*v_deltacd\n",
      "variable deltacdm2A1 equal -0.115*v_deltacd\n",
      "\n",
      "thermo ${freq\n",
      "\n",
      "# Set-up Loo\n",
      "variable runid loop 1 ${nblocks} pa\n",
      "    label runloop\n",
      "\n",
      "    # Adjust param for the box and equilibrat\n",
      "    variable param equal v_lambdas[v_runid\n",
      "    if \"${runid} == 1\" then \n",
      "        \"jump SELF skipequil\n",
      "    variable ind equal v_runid-\n",
      "    variable param0 equal v_lambdas[v_ind\n",
      "    variable paramramp equal ramp(v_param0,v_param\n",
      "    variable paramrampA0 equal -0.115*ramp(v_param0,v_param\n",
      "    variable paramA0 equal -0.115*v_para\n",
      "    variable paramrampA1 equal 0.115*ramp(v_param0,v_param\n",
      "    variable paramA1 equal 0.115*v_para\n",
      "    fix ADAPT all adapt/fep 1 \n",
      "        atom charge 1  v_paramrampA0 \n",
      "        atom charge 2  v_paramrampA1\n",
      "    thermo_style custom step v_paramramp v_paramrampA0 v_paramrampA1 temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Ram\n",
      "    unfix ADAP\n",
      "    fix ADAPT all adapt/fep ${freq} \n",
      "        atom charge 1  v_paramA0 \n",
      "        atom charge 2  v_paramA1\n",
      "    thermo_style custom step v_param v_paramA0 v_paramA1 temp press pe evdwl enthalp\n",
      "    run ${runtime_equil} # Run Equi\n",
      "\n",
      "    label skipequil\n",
      "    write_data files/charge_${param}.dat\n",
      "\n",
      "    # Initialize compute\n",
      "    thermo_style custom step v_param temp press pe evdwl enthalp\n",
      "    compute FEPdb all fep ${TK} \n",
      "        atom charge 1  v_deltacdm2A0 \n",
      "        atom charge 2  v_deltacdm2A1 \n",
      "        volume ye\n",
      "    compute FEPdf all fep ${TK} \n",
      "        atom charge 1  v_deltacdmA0 \n",
      "        atom charge 2  v_deltacdmA1 \n",
      "        volume ye\n",
      "    fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_tinst v_pinst v_vinst v_pe \n",
      "        c_FEPdb[1] c_FEPdf[1] file files/ti_charge_${param}.txt\n",
      "    variable delta0 equal v_lambdas[1]-v_para\n",
      "    variable deltaA00 equal -0.115*v_delta\n",
      "    variable deltaA01 equal 0.115*v_delta\n",
      "    compute FEP000 all fep ${TK} \n",
      "        atom charge 1  v_deltaA00 \n",
      "        atom charge 2  v_deltaA01 \n",
      "        volume ye\n",
      "    variable param000 equal v_param+v_delta\n",
      "    fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_pe \n",
      "        c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_charge_${param}_${param000}.txt\n",
      "    variable delta1 equal v_lambdas[2]-v_para\n",
      "    variable deltaA10 equal -0.115*v_delta\n",
      "    variable deltaA11 equal 0.115*v_delta\n",
      "    compute FEP001 all fep ${TK} \n",
      "        atom charge 1  v_deltaA10 \n",
      "        atom charge 2  v_deltaA11 \n",
      "        volume ye\n",
      "    variable param001 equal v_param+v_delta\n",
      "    fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_pe \n",
      "        c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_charge_${param}_${param001}.txt\n",
      "    variable delta2 equal v_lambdas[3]-v_para\n",
      "    variable deltaA20 equal -0.115*v_delta\n",
      "    variable deltaA21 equal 0.115*v_delta\n",
      "    compute FEP002 all fep ${TK} \n",
      "        atom charge 1  v_deltaA20 \n",
      "        atom charge 2  v_deltaA21 \n",
      "        volume ye\n",
      "    variable param002 equal v_param+v_delta\n",
      "    fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_pe \n",
      "        c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_charge_${param}_${param002}.txt\n",
      "    variable delta3 equal v_lambdas[4]-v_para\n",
      "    variable deltaA30 equal -0.115*v_delta\n",
      "    variable deltaA31 equal 0.115*v_delta\n",
      "    compute FEP003 all fep ${TK} \n",
      "        atom charge 1  v_deltaA30 \n",
      "        atom charge 2  v_deltaA31 \n",
      "        volume ye\n",
      "    variable param003 equal v_param+v_delta\n",
      "    fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_pe \n",
      "        c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_charge_${param}_${param003}.txt\n",
      "    variable delta4 equal v_lambdas[5]-v_para\n",
      "    variable deltaA40 equal -0.115*v_delta\n",
      "    variable deltaA41 equal 0.115*v_delta\n",
      "    compute FEP004 all fep ${TK} \n",
      "        atom charge 1  v_deltaA40 \n",
      "        atom charge 2  v_deltaA41 \n",
      "        volume ye\n",
      "    variable param004 equal v_param+v_delta\n",
      "    fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_pe \n",
      "        c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_charge_${param}_${param004}.txt\n",
      "    variable delta5 equal v_lambdas[6]-v_para\n",
      "    variable deltaA50 equal -0.115*v_delta\n",
      "    variable deltaA51 equal 0.115*v_delta\n",
      "    compute FEP005 all fep ${TK} \n",
      "        atom charge 1  v_deltaA50 \n",
      "        atom charge 2  v_deltaA51 \n",
      "        volume ye\n",
      "    variable param005 equal v_param+v_delta\n",
      "    fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_pe \n",
      "        c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_charge_${param}_${param005}.txt\n",
      "  #  dump TRAJ all custom ${freq} files/dump_charge_${param}.lammpstrj id mol type element xu yu z\n",
      "\n",
      "    run ${runtime_prod}\n",
      "    uncompute FEPd\n",
      "    uncompute FEPd\n",
      "    if \"${runid} != 1\" then \n",
      "        \"unfix ADAPT\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    uncompute FEP00\n",
      "    unfix FEPout00\n",
      "    unfix FEPou\n",
      "  #  undump TRA\n",
      "\n",
      "    next runi\n",
      "    jump SELF runloop\n",
      "thermo_style custom v_param temp press pe evdwl enthalp\n",
      "write_data final.data nocoef\n"
     ]
    }
   ],
   "source": [
    "# Input where charges change\n",
    "output = gali.generate_mbar_input(\n",
    "    [0.0, 1.0], # Scale the charges from zero to in full force\n",
    "    0.2, \n",
    "    [\n",
    "        [\"atom\", (\"charge\", 1, -0.115000)],\n",
    "        [\"atom\", (\"charge\", 2, 0.115000)],\n",
    "    ],\n",
    "    300, \n",
    "#    output_file=\"your/output/here\", # Optional file to write to\n",
    "    n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)\n",
    "    n_run_prod_steps=5000000,  # Duration of production run for each step (i.e., window)\n",
    ")\n",
    "for line in output:\n",
    "    print(line[:-2])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gali",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}