################ DESCRIPTION ###############

Example for how to generate an run a simulation with the All-atom model using SMOG 2 and OpenSMOG
In this example, we will add a coulomb potential between charges, and we will add a +1 charge to each N in LYS. We will further add some K+, Cl- and MG2+ ions. We will also add some lines to the templates, so that they can be used with OpenSMOG.


############### PREREQUISITES ##############

- SMOG 2 must be configured and all SMOG tools must be in your PATH.  If you are using a SMOG container, this will be satisfied.
- the alias "python" must map to a version for which the OpenSMOG libraries are available. Conda environments can be helpful to set this up. 

 
############## PRE-SMOG2 STEPS #############

Copy the default AA model templates from <SMOG 2 DIR>/share/templates/SBM_AA
To avoid confusion, rename all files SBM_AA+coulomb+ions
In the .bif file, find the residue entry for LYS. 

  Under LYS, replace the following line:
<atom bType="B_1" nbType="NB_1" pairType="C_1">NZ</atom>

  with
<atom bType="B_1" nbType="NB_1" pairType="C_1" charge="1" >NZ</atom>


In the .sif file, under the <sif> element, list the following child element:
   <OpenSMOGsettings>
    <constants>
      <!-- Note: K_coul is 1/5 of the value in most software
                 since SMOG models use reduced units. For discussion
                 on units in this model, see SMOG-ion manuscript
      -->
      <constant name="K_coul" value="27.787097"/>
      <constant name="dielectric" value="80"/>
      <!-- Define all atoms to be the same size. This is necessary
                 since OpenSMOG nonbonded parameters overwrite
                 anything that may have been found in the .top file--> 
      <constant name="atomsigma" value="0.25"/>
      <constant name="atomepsilon" value="0.1"/>
    </constants>
  </OpenSMOGsettings>

and 
  <CustomNonBonded 
        OpenSMOGparameters=" " 
        OpenSMOGcombrule=" none " 
        OpenSMOGpotential="K_coul*q1*q2/dielectric*(1/r-1/r_c)+atomepsilon*atomsigma^12*(1/r^12-1/r_c^12)"
  />

Explainer: For convenience, we defined constants that could be used in our force field definition.  We then defined a generic Coulomb+1/r^12 term (truncated) for all nonbonded interactions. This is not the only way one could go about introducing this potential, but it will work. It also places all parameters in a single location, as opposed to having to specify each at various stages of simulation preparation.

Finally, we will also create a file called SBM_AA+coulomb+ions.ions.def, which will have the contents:
; name  mass   charge   C12             C6 
MG      1       2       5.96046e-09     0
K       1       1       5.96046e-09     0
CL      1       -1      5.96046e-09     0


######### USING SMOG2 AND SMOG-TOOLS #######
 
Begin with the file 2ci2.pdb, obtained directly from RCSB. For your convenience, it is in the current directory.

Preprocess the PDB file: smog_adjustPDB will remove water molecules and recognize all valid PDB keywords. Accordingly, smog_adjustPDB may be applied directly to the downloaded PDB file.
Since this example set of templates uses naming conventions that are consistent with the default SMOG model, we do not need to specify the mapping file.  The default mapping rules will be applied.
 
>smog_adjustPDB -i 2ci2.pdb -o 2ci2.adjusted.pdb -removewater



# Generate your force field for the biomolecule with smog2:  
- The -t flag indicates that the template directory SBM_AA+coulomb+ions should be used to define the force field.
- The -OpenSMOG flag indicates that the input files should be prepared for use with OpenSMOG/OpenMM
- The option "-boxbuffer 5" will make a box that extends beyond the protein by 5nm in each direction
- The -center option will center the protein in the box (equidistant from boundaries)

>smog2 \
	-i 2ci2.adjusted.pdb \
	-t SBM_AA+coulomb+ions \
	-OpenSMOG \
	-dname 2ci2.OpenSMOG.AA+coulomb+ions \
	-center \
	-boxbuffer 5

For this system, you should receive a message indicating that the total charge is +6.

Now we need to add some ions. The ion definitions are given in the template file SBM_AA+coulomb+ions.ions.def.  The file gives the ion name, mass, charge, c12 and c6 parameters.

To add ions, we will use the smog_ions tool. Here, we are going to define the ions based on the *.ions.def file in the template directory. We would like to add 2 MG2+ ions

>smog_ions \
	-t SBM_AA+coulomb+ions \
	-f 2ci2.OpenSMOG.AA+coulomb+ions.top \
	-g 2ci2.OpenSMOG.AA+coulomb+ions.gro \
	-OpenSMOG 2ci2.OpenSMOG.AA+coulomb+ions.xml \
	-of 2ci2.OpenSMOG.AA+coulomb+ions.MG.top \
	-og 2ci2.OpenSMOG.AA+coulomb+ions.MG.gro \
	-OpenSMOGout 2ci2.OpenSMOG.AA+coulomb+ions.MG.xml \
	-ionnm MG 
	-ionn 2

By adding 2 Mg2+, we bring the total charge to 10.

We will now neutralize the system with 10 Cl- ions.

>smog_ions \
	-t SBM_AA+coulomb+ions \
	-f 2ci2.OpenSMOG.AA+coulomb+ions.MG.top \
	-g 2ci2.OpenSMOG.AA+coulomb+ions.MG.gro \
	-OpenSMOG 2ci2.OpenSMOG.AA+coulomb+ions.MG.xml \
	-of 2ci2.OpenSMOG.AA+coulomb+ions.MGCL.top \
	-og 2ci2.OpenSMOG.AA+coulomb+ions.MGCL.gro \
	-OpenSMOGout 2ci2.OpenSMOG.AA+coulomb+ions.MGCL.xml  \
	-ionnm CL \
	-ionn 10


#### RUNNING A SIMULATION WITH OPENSMOG ####

# At this point, you have all files needed to start your calculations in OpenSMOG.  Using the script run.ci2.OpenSMOG.AA+coulomb+ions.py (in this directory), you can load, energy minimize and simulate the system. See comments in .py file for explanation of settings.

>python run.ci2.OpenSMOG.AA+coulomb+ions.py