.. _prepare:

How to Prepare Your Data
========================

Preparing your biomolecular system for Wepy simulation mostly follows
the same steps you would use to run a normal molecular dynamics (MD)
trajectory with OpenMM.

To prepare your data for Wepy you can follow these 3 steps:

#. Use a system builder (e.g. CharmmGUI) to build your system.
#. Perform a short MD simulation with OpenMM.
#. Use the simulation output to build a restart file, as well as pickle
   files for the system and topology.

We describe a protocol using CharmmGUI below, but this can be easily
adapted to accomodate other workflows as well.

Building your system with CharmmGUI
-----------------------------------

`CharmmGUI <https://www.charmm-gui.org/>`__ is a powerful web-based
interface for preparing input files for MD simulations. There are a
number of resources on that page to assist you in choosing the
appropriate options for your project. A minimal workflow is provided
below:

#. **Visit CharmmGUI:** Open your web browser and go to
   https://www.charmm-gui.org . You will need to create an account
   and/or log in.

#. **Select the Appropriate Module:** Choose the module that best suits
   your needs (e.g., "Membrane Builder", "Solution Builder", etc.). For
   most MD simulations, the "Solution Builder" is commonly used.

#. **Upload Your Molecular Structure**: Upload your molecular structure
   file (e.g., PDB file) to the CharmmGUI interface, or enter the PDB
   ID. If uploading, ensure your file is correctly formatted and
   contains all necessary components.

#. **Set Up Your Simulation Parameters:** Follow the guided steps in
   CharmmGUI to set up the details of your simulation. This typically
   includes selecting the appropriate chains and entities, defining the
   type of solvent, adding ions, setting the box size, temperature and
   other relevant details.

#. **Generate Input Files:** Once you have configured your system,
   CharmmGUI will generate the necessary input files for your MD
   simulation. **For output type, be sure to select OpenMM**. This will
   include all of the necessary topology files, coordinate files, and
   parameter files.

#. **Download the Input Files:** Download the tgz-format archive file to
   your local machine and/or move them into your cluster or computing
   environment. Expand this using:

   ``tar -xzf <charmm-gui-file>.tgz``.

   The files are now in the directory ``charmm-gui-XXXXXXXXXX/openmm``
   where XXXXXXXXXX is your particular job number. You will use the
   files in this directory to perform the MD simulation in the next
   step.

MD Simulation
-------------

The output from CharmmGUI is not yet equilibrated. It should be energy
minimized and then gradually heated to the production temperature.
CharmmGUI has provided scripts for this purpose, which are contained in
the ``README`` file.

A nice way to carry out this equilibration is to make this file
executable (with ``chmod +x README``) and then execute its contents:

.. code:: bash

   ./README

This can also be run in a SLURM script, if working on a compute cluster:

.. code:: bash

   #!/bin/bash
   #SBATCH --job-name=charmm-gui-equil
   #SBATCH --gpus=1
   #SBATCH --time=4:00:00
   #SBATCH --mem=5G

   # activate the conda environment with openmm
   conda activate wepy-env

   # load a CUDA module if necessary
   module load CUDA/11.8.0

   cd ~/charmm-gui-XXXXXXXXXX/openmm

   ./README >& t.log 

Building System and Topology Pickle Files
-----------------------------------------

Now that you have your equilibrated system, you can use the simulation
output to build system and topology pickle files. To achieve this you
will use the ``openmm_run.py`` script in the
``charmm-gui-XXXXXXXXXX/openmm`` folder and the corresponding README
file generated by CHARMM-GUI.

To build the system and topology pickle files, follow these steps:

#. **Modify the openmm\ run.py script:** This script is generated by
   CHARMM-GUI. It contains the necessary code to generate the system and
   topology objects which we will save in pkl files, and then exit the
   script. An example of this modified file is given below, which is
   also included in ``examples/charmm-gui/openmm_make_pkls.py``.

   .. code:: python

      """
      Generated by CHARMM-GUI (http://www.charmm-gui.org)

      openmm_run.py

      This program is OpenMM running scripts written in python.

      Correspondance: jul316@lehigh.edu or wonpil@lehigh.edu
      Last update: June 18, 2021
      """

      from __future__ import print_function
      import argparse
      import sys
      import os
      import pickle as pkl

      from omm_readinputs import *
      from omm_readparams import *
      from omm_vfswitch import *
      from omm_barostat import *
      from omm_restraints import *
      from omm_rewrap import *

      from openmm.unit import *
      from openmm import *
      from openmm.app import *

      parser = argparse.ArgumentParser()
      parser.add_argument('--platform', nargs=1, help='OpenMM platform (default: CUDA or OpenCL)')
      parser.add_argument('-i', dest='inpfile', help='Input parameter file', required=True)
      parser.add_argument('-p', dest='topfile', help='Input topology file', required=True)
      parser.add_argument('-c', dest='crdfile', help='Input coordinate file', required=True)
      parser.add_argument('-t', dest='toppar', help='Input CHARMM-GUI toppar stream file (optional)')
      parser.add_argument('-b', dest='sysinfo', help='Input CHARMM-GUI sysinfo stream file (optional)')
      parser.add_argument('-ff', dest='fftype', help='Input force field type (default: CHARMM)', default='CHARMM')
      parser.add_argument('-icrst', metavar='RSTFILE', dest='icrst', help='Input CHARMM RST file (optional)')
      parser.add_argument('-irst', metavar='RSTFILE', dest='irst', help='Input restart file (optional)')
      parser.add_argument('-ichk', metavar='CHKFILE', dest='ichk', help='Input checkpoint file (optional)')
      parser.add_argument('-opdb', metavar='PDBFILE', dest='opdb', help='Output PDB file (optional)')
      parser.add_argument('-orst', metavar='RSTFILE', dest='orst', help='Output restart file (optional)')
      parser.add_argument('-ochk', metavar='CHKFILE', dest='ochk', help='Output checkpoint file (optional)')
      parser.add_argument('-odcd', metavar='DCDFILE', dest='odcd', help='Output trajectory file (optional)')
      parser.add_argument('-rewrap', dest='rewrap', help='Re-wrap the coordinates in a molecular basis (optional)', action='store_true', default=False)
      args = parser.parse_args()

      # Load parameters
      print("Loading parameters")
      inputs = read_inputs(args.inpfile)

      top = read_top(args.topfile, args.fftype.upper())
      crd = read_crd(args.crdfile, args.fftype.upper())
      if args.fftype.upper() == 'CHARMM':
      params = read_params(args.toppar)
      top = read_box(top, args.sysinfo) if args.sysinfo else gen_box(top, crd)

      # Build system
      nboptions = dict(nonbondedMethod=inputs.coulomb,
                     nonbondedCutoff=inputs.r_off*nanometers,
                     constraints=inputs.cons,
                     ewaldErrorTolerance=inputs.ewald_Tol)
      if inputs.vdw == 'Switch': nboptions['switchDistance'] = inputs.r_on*nanometers
      if inputs.vdw == 'LJPME':  nboptions['nonbondedMethod'] = LJPME
      if inputs.implicitSolvent:
      nboptions['implicitSolvent'] = inputs.implicitSolvent
      nboptions['implicitSolventSaltConc'] = inputs.implicit_salt*(moles/liter)
      nboptions['temperature'] = inputs.temp*kelvin
      nboptions['soluteDielectric'] = inputs.solut_diele
      nboptions['solventDielectric'] = inputs.solve_diele
      nboptions['gbsaModel'] = inputs.gbsamodel

      if   args.fftype.upper() == 'CHARMM': system = top.createSystem(params, **nboptions)
      elif args.fftype.upper() == 'AMBER':  system = top.createSystem(**nboptions)

      if inputs.vdw == 'Force-switch': system = vfswitch(system, top, inputs)
      if inputs.lj_lrc == 'yes':
      for force in system.getForces():
            if isinstance(force, NonbondedForce): force.setUseDispersionCorrection(True)
               if isinstance(force, CustomNonbondedForce) and force.getNumTabulatedFunctions() != 1:
                  force.setUseLongRangeCorrection(True)
      if inputs.e14scale != 1.0:
      for force in system.getForces():
            if isinstance(force, NonbondedForce): nonbonded = force; break
      for i in range(nonbonded.getNumExceptions()):
               atom1, atom2, chg, sig, eps = nonbonded.getExceptionParameters(i)
               nonbonded.setExceptionParameters(i, atom1, atom2, chg*inputs.e14scale, sig, eps)

      if inputs.pcouple == 'yes':      system = barostat(system, inputs)
      if inputs.rest == 'yes':         system = restraints(system, crd, inputs)

      # system and topology are made!
      # write the pkl files
      pkl.dump(system,open('system.pkl','wb'))
      pkl.dump(top.topology,open('topology.pkl','wb'))

   NOTE again that this is just a truncated version of the
   ``openmm_run.py`` script with two ``pkl.dump`` commands at the end.
   It will create system and topology objects that correspond to the
   options specified in the CharmmGUI configuration files.

#. **Call the new script as done in the README:** The README file, also
   generated by CHARMM-GUI, provides the call signature for running this
   script with the "production" parameters. Check your README script for
   the appropriate call signature. A sample command is as follows:

   .. code:: bash

      python -u openmm_make_pkls.py -i step5_production.inp -t toppar.str -p step3_input.psf -c step3_input.crd 

   This script will generate two pickle files: ``system.pkl`` and
   ``topology.pkl``. These files will be used in the next step to run
   the simulation with the generated system and topology. Now you are
   all set to start setting up the Wepy simulation.

​*\* For rebinding simulations, you will need to generate system and
topology files for the bonded state as well. \*\*