<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Proton Titration of a Rigid Body in Explicit Salt Solution
This will simulate proton equilibria in a rigid body composed of particles in a spherical simulation container. We use a continuum solvent and explicit, soft spheres for salt particles which are treated grand canonically. During simulation, titratable sites are updated with swap moves according to input pH.
System Requirements
This Jupyter Notebook was originally run in MacOS 10.11 with GCC 4.8, Python2, matplotlib, pandas within the Anaconda environment. Contemporary Linux distributions such as Ubuntu 14.04 should work as well.
Download and build Faunus
We use a custom Metropolis Monte Carlo (MC) program build within the Faunus framework. The sections below will fetch the complete faunus project and compile the program.
Step1:  Let's coarse grain an atomic PDB structure to the amino acid level
For this purpose we use the module pytraj and also locate possible disulfide bonds as these should not be allowed to titrate. Note that if found, the residues in the generated .aam file should be manually renamed to prevent them from being regarded as titratable. Likewise, the N and C terminal ends need to be handled manually.
Step2:  Create Input and run MC simulation
Step3:  Analysis
<ASSISTANT_TASK:>
Python Code:
from __future__ import division, unicode_literals, print_function
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np, pandas as pd
import os.path, os, sys, json, filecmp, copy
plt.rcParams.update({'font.size': 16, 'figure.figsize': [8.0, 6.0]})
try:
    workdir
except NameError:
    workdir=%pwd
else:
    %cd $workdir
print(workdir)
%%bash -s "$workdir"
%cd $1
if [ ! -d "faunus/" ]; then
    git clone https://github.com/mlund/faunus.git
    cd faunus
    git checkout 86a1f74
else
    cd faunus
fi
# if different, copy custom gctit.cpp into faunus
if ! cmp ../titrate.cpp src/examples/gctit.cpp
then
    cp ../titrate.cpp src/examples/gctit.cpp
fi
pwd
cmake . -DCMAKE_BUILD_TYPE=Release -DENABLE_APPROXMATH=on &>/dev/null
make example_gctit -j4
Explanation: Proton Titration of a Rigid Body in Explicit Salt Solution
This will simulate proton equilibria in a rigid body composed of particles in a spherical simulation container. We use a continuum solvent and explicit, soft spheres for salt particles which are treated grand canonically. During simulation, titratable sites are updated with swap moves according to input pH.
System Requirements
This Jupyter Notebook was originally run in MacOS 10.11 with GCC 4.8, Python2, matplotlib, pandas within the Anaconda environment. Contemporary Linux distributions such as Ubuntu 14.04 should work as well.
Download and build Faunus
We use a custom Metropolis Monte Carlo (MC) program build within the Faunus framework. The sections below will fetch the complete faunus project and compile the program.
End of explanation
%cd $workdir
import mdtraj as md
traj = md.load_pdb('1BPI.pdb')
for chain in traj.topology.chains:
    print('chain: ', chain.index)
    
    # filter pdb to only select protein(s)
    sel = chain.topology.select('protein')
    top = chain.topology.subset(sel)
    f = open('chain'+str(chain.index)+'.aam','w')
    f.write(str(top.n_residues)+'\n')
    
    # locate disulfide bonds (not used for anything yet)
    for b in top.bonds:
        i = b[0].residue.index
        j = b[1].residue.index
        if j>i+1:
            if (b[0].residue.name == 'CYS'):
                if (b[1].residue.name == 'CYS'):
                    print('SS-bond between residues', i, j)
        
    # loop over residues and calculate residue mass centers, radius, and weight
    top.create_disulfide_bonds( traj.xyz[0] )
    for res in top.residues:
        if res.is_protein:
            cm = [0,0,0] # residue mass center
            mw = 0       # residue weight
            for a in res.atoms:
                cm = cm + a.element.mass * traj.xyz[0][a.index]
                mw = mw + a.element.mass
            cm = cm/mw*10
            radius = ( 3./(4*np.pi)*mw/1.0 )**(1/3.)
            f.write('{0:4} {1:5} {2:8.3f} {3:8.3f} {4:8.3f} {5:6.3f} {6:6.2f} {7:6.2f}\n'\
                    .format(res.name,res.index,cm[0],cm[1],cm[2],0,mw,radius))
    f.close()
Explanation: Let's coarse grain an atomic PDB structure to the amino acid level
For this purpose we use the module pytraj and also locate possible disulfide bonds as these should not be allowed to titrate. Note that if found, the residues in the generated .aam file should be manually renamed to prevent them from being regarded as titratable. Likewise, the N and C terminal ends need to be handled manually.
End of explanation
pH_range       = [7.0]
salt_range     = [0.03] # mol/l
%cd $workdir'/'
def mkinput():
    js = {
        "energy": {
            "eqstate": { "processfile": "titrate.json" },
            "nonbonded": {
                "coulomb": { "epsr": 80 }
            }
        },
        "system": {
            "temperature": 298.15,
            "sphere" : { "radius" : 90 },
            "mcloop": { "macro": 10, "micro": micro }
        },
        "moves": {
            "gctit"         : { "molecule": "salt", "prob": 0.5 },
            "atomtranslate" : {
                "salt":  { "prob": 0.5 }
            }
        },
        "moleculelist": {
            "protein":  { "structure":"../chain0.aam", "Ninit":1, "insdir":"0 0 0"},
            "salt": {"atoms":"Na Cl", "Ninit":60, "atomic":True }
        },
        "atomlist" : {
            "Na"   :  { "q": 1, "r":1.9, "eps":0.005, "mw":22.99, "dp":100, "activity":salt },
            "Cl"   :  { "q":-1, "r":1.7, "eps":0.005, "mw":35.45, "dp":100, "activity":salt },
            "ASP"  :  { "q":-1, "r":3.6, "eps":0.05, "mw":110 },
            "HASP" :  { "q":0,  "r":3.6, "eps":0.05, "mw":110 },
            "LASP" :  { "q":2,  "r":3.6, "eps":0.05, "mw":110 },
            "CTR"  :  { "q":-1, "r":2.0, "eps":0.05, "mw":16 },
            "HCTR" :  { "q":0,  "r":2.0, "eps":0.05, "mw":16 },
            "GLU"  :  { "q":-1, "r":3.8, "eps":0.05, "mw":122 },
            "HGLU" :  { "q":0,  "r":3.8, "eps":0.05, "mw":122 },
            "LGLU" :  { "q":2,  "r":3.8, "eps":0.05, "mw":122 },
            "HIS"  :  { "q":0,  "r":3.9, "eps":0.05, "mw":130 },
            "HHIS" :  { "q":1,  "r":3.9, "eps":0.05, "mw":130 },
            "NTR"  :  { "q":0,  "r":2.0, "eps":0.05, "mw":14 },
            "HNTR" :  { "q":1,  "r":2.0, "eps":0.05, "mw":14 },
            "TYR"  :  { "q":-1, "r":4.1, "eps":0.05, "mw":154 },
            "HTYR" :  { "q":0,  "r":4.1, "eps":0.05, "mw":154 },
            "LYS"  :  { "q":0,  "r":3.7, "eps":0.05, "mw":116 },
            "HLYS" :  { "q":1,  "r":3.7, "eps":0.05, "mw":116 },
            "CYb"  :  { "q":0,  "r":3.6, "eps":0.05, "mw":103 },
            "CYS"  :  { "q":-1, "r":3.6, "eps":0.05, "mw":103 },
            "HCYS" :  { "q":0,  "r":3.6, "eps":0.05, "mw":103 },
            "ARG"  :  { "q":0,  "r":4.0, "eps":0.05, "mw":144 },
            "HARG" :  { "q":1,  "r":4.0, "eps":0.05, "mw":144 },
            "ALA"  :  { "q":0,  "r":3.1, "eps":0.05, "mw":66 },
            "ILE"  :  { "q":0,  "r":3.6, "eps":0.05, "mw":102 },
            "LEU"  :  { "q":0,  "r":3.6, "eps":0.05, "mw":102 },
            "MET"  :  { "q":0,  "r":3.8, "eps":0.05, "mw":122 },
            "PHE"  :  { "q":0,  "r":3.9, "eps":0.05, "mw":138 },
            "PRO"  :  { "q":0,  "r":3.4, "eps":0.05, "mw":90 },
            "TRP"  :  { "q":0,  "r":4.3, "eps":0.05, "mw":176 },
            "VAL"  :  { "q":0,  "r":3.4, "eps":0.05, "mw":90 },
            "SER"  :  { "q":0,  "r":3.3, "eps":0.05, "mw":82 },
            "THR"  :  { "q":0,  "r":3.5, "eps":0.05, "mw":94 },
            "ASN"  :  { "q":0,  "r":3.6, "eps":0.05, "mw":108 },
            "GLN"  :  { "q":0,  "r":3.8, "eps":0.05, "mw":120 },
            "GLY"  :  { "q":0,  "r":2.9, "eps":0.05, "mw":54 }
        },
        "processes" : {
            "H-Asp" : { "bound":"HASP" , "free":"ASP" , "pKd":4.0  , "pX":pH },
            "H-Ctr" : { "bound":"HCTR" , "free":"CTR" , "pKd":2.6  , "pX":pH },
            "H-Glu" : { "bound":"HGLU" , "free":"GLU" , "pKd":4.4  , "pX":pH },
            "H-His" : { "bound":"HHIS" , "free":"HIS" , "pKd":6.3  , "pX":pH },
            "H-Arg" : { "bound":"HARG" , "free":"ARG" , "pKd":12.0 , "pX":pH },
            "H-Ntr" : { "bound":"HNTR" , "free":"NTR" , "pKd":7.5  , "pX":pH },
            "H-Cys" : { "bound":"HCYS" , "free":"CYS" , "pKd":10.8 , "pX":pH },
            "H-Tyr" : { "bound":"HTYR" , "free":"TYR" , "pKd":9.6  , "pX":pH },
            "H-Lys" : { "bound":"HLYS" , "free":"LYS" , "pKd":10.4 , "pX":pH }
        }
    }
    with open('titrate.json', 'w+') as f:
        f.write(json.dumps(js, indent=4))
for pH in pH_range:
    for salt in salt_range:
        pfx='pH'+str(pH)+'-I'+str(salt)
        if not os.path.isdir(pfx):
            %mkdir -p $pfx 
            %cd $pfx
            # equilibration run (no translation)
            !rm -fR state
            micro=100
            mkinput()
            !../faunus/src/examples/gctit > eq
            # production run
            micro=1000
            mkinput()
            %time !../faunus/src/examples/gctit > out
            
            
            %cd ..
%cd ..
print('done.')
Explanation: Create Input and run MC simulation
End of explanation
%cd $workdir'/'
import json
for pH in pH_range:
    for salt in salt_range:
        pfx='pH'+str(pH)+'-I'+str(salt)
        if os.path.isdir(pfx):
            %cd $pfx            
            js = json.load( open('gctit-output.json') )
            charge = js['protein']['charge']
            index = js['protein']['index']
            resname = js['protein']['resname']
            plt.plot(index,charge, 'ro')           
            %cd ..
for i in range(0,len(index)):
    label = resname[i]+' '+str(index[i]+1)
    plt.annotate(label, xy=(index[i], charge[i]), fontsize=8, rotation=70)
plt.title('Protonation States of All Residues')
plt.legend(loc=0, frameon=False)
plt.xlabel(r'residue number')
plt.ylabel(r'average charge, $z$')
plt.ylim((-1.1, 1.1))
#plt.xticks(i, resname, rotation=70, fontsize='small')
plt.savefig('fig.pdf', bbox_inches='tight')
Explanation: Analysis
End of explanation
<END_TASK>