<jupyter_start><jupyter_text>Worksheet 8: Coarse-grained modeles<jupyter_code># put import statements here
import numpy as np
import mdtraj as md
import simtk.openmm as mm

# Force field parametrization constants
ATOM_MASS = 1.0         # a.u.
BOND_STRENGTH = 20000.0 # kJ/(nm**2 mol)
ANGLE_STRENGTH = 40.0   # kJ/(mol rad**2)
NONNATIVE_SIGMA = 0.4   # nm, excluded radius for nonnative interactions

# CG modeling parameters
TEMPERATURE = 132   # K,  modeling temperature
STRIDE = 1000       # steps, record coordinates every STRIDE steps
STEPSIZE = 0.0005   # ps,  integration step
FRICTION = 1        # 1/ps, friction coefficient for LangevinIntegrator
# Input files
pdb_code = '1PGB'
aa_pdb_file = f'{pdb_code}.pdb'  # aa stands for all-atom
native_contacts_file = '1PGB_native_contacts.dat' # Each row contains the pair of native contacts, indexed from 0<jupyter_output><empty_output><jupyter_text>Task 1: Model building  and simulation<jupyter_code>native_contacts  = np.loadtxt(native_contacts_file,dtype=int)
def add_nonbonded_interactions(ca_structure, system, native_contacts):
    """
    Add nonbonded interactions to OpenMM System object
    
    Parameters
    ----------
    ca_structure :  mdtraj.Trajectory object
        Trajectory that contains native structure for the CG model.
        If more than one frame provided, only the first frame will be used
        and the rest will be discarded. 
    
    native_contacts : numpy ndarray, dtype int
        A Mx2 numpy array. Each row  represents a pair of residues that are in contact in the 
        native structure
        
    system : simtk.openmm.openmm.System object
        Object to which forces should be added
        
    Returns
    -------
    system : simtk.openmm.openmm.System object
       Modified system
    
    """
    # Add nonbonded interaction 
    n_atoms = ca_structure.top.n_atoms
    nonbonded_force = mm.openmm.CustomNonbondedForce("(sigma/r)^12; sigma=sqrt(sigma1*sigma2)")
    nonbonded_force.addPerParticleParameter("sigma")
    for i in range(n_atoms):
        nonbonded_force.addParticle([NONNATIVE_SIGMA])
        
    # Exclude 1-2, 1-3 and 1-4 interactions:
    # Bond cutoff=3: atoms separated by 3 bonds or less are excluded
    nonbonded_force.createExclusionsFromBonds(bonds=[[i, i+1] for i in range(n_atoms-1)], bondCutoff=3)
            
    # Add native interaction
    native_force = mm.openmm.CustomBondForce("5.0*((rnot/r)^12) - 6.0*((rnot/r)^10)")
    native_force.addPerBondParameter("rnot")
    distances =  md.compute_distances(ca_structure, native_contacts)[0]
    for ndx, contact in enumerate(native_contacts):
        # OpenMM does not work correctly with numpy objects, all scalars should be 
        # explicitly converted
        atom_1 = int(contact[0])
        atom_2 = int(contact[1])
        native_force.addBond(atom_1, atom_2, [float(distances[ndx])])
        nonbonded_force.addExclusion(atom_1, atom_2)
        
    # Add forces to the system
    system.addForce(nonbonded_force)
    system.addForce(native_force)
    return system 
                                     
def add_dihedral_interactions(ca_structure, system):
    """
    Add dihedral interactions to OpenMM System object
    
    Parameters
    ----------
    ca_structure :  mdtraj.Trajectory object
        Trajectory that contains native structure for the CG model.
        If more than one frame provided, only the first frame will be used
        and the rest will be discarded. 
        
    system : simtk.openmm.openmm.System object
        Object to which forces should be added
        
    Returns
    -------
    system : simtk.openmm.openmm.System object
       Modified system

    """
    # Add dihedral angle potential. 
    atom_quadruples = [[i, i+1, i+2, i+3] for i in range(n_atoms - 3)]
    dihedrals = md.compute_dihedrals(ca_structure, atom_quadruples)[0] 
    dihedral_force = mm.openmm.PeriodicTorsionForce()
    for dihedral_ndx in range(n_atoms - 3):
        dihedral_1 = dihedrals[dihedral_ndx] + np.pi
        dihedral_2 = 3*(dihedrals[dihedral_ndx] + np.pi)
        dihedral_force.addTorsion(*atom_quadruples[dihedral_ndx], 1, dihedral_1, 1.0)
        dihedral_force.addTorsion(*atom_quadruples[dihedral_ndx], 3, dihedral_2, 0.5)
    system.addForce(dihedral_force)
    return system


def save_system(system, name='OpenMM_system.xml'):
    """
    Save OpenMM system object to a human readable xml file
    
    Parameters:
    ----------- 
    
    name : str, optional
        Name of the file, should end with .xml extension
        
    system : simtk.openmm.openmm.System object
    """
    with open(name, 'w') as f:
        system_xml = mm.XmlSerializer.serialize(system)
        f.write(system_xml) 
    return
# <add your code here><jupyter_output><empty_output><jupyter_text>Task 2 Analysis<jupyter_code>def get_fraction_native_contacts(traj, native_contacts, cutoff=1.15):
    """
    Calculate fraction of formed native contacts for each frame in trajectory
    
    Parameters
    ----------
     traj :  mdtraj.Trajectory object
     
    native_contacts : numpy ndarray, dtype int
        A Mx2 numpy array. Each row  represents a pair of residues that are in contact in the 
        native structure
        
    cutoff : float, optional
        If two beads are separated by distance lower than `cutoff`, the contact is considered formed.
        
    Returns
    -------
    
    fraction_of_contacts : 1d numpy ndarray
        element i of the array represents fraction of native contact formed in that frame i of `traj`.

    """
    distances = md.compute_distances(traj, native_contacts)
    contact_formed = np.where(distances < cutoff, 1, 0)
    num_of_contacts = np.sum(contact_formed, axis=1)
    fraction_of_contacts = num_of_contacts/contact_formed.shape[1]
    return(fraction_of_contacts)
# <add your code here><jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0855497204811113, 'fraction_numerical': 0.015924106386583092, 'mean_word_length': 2.92814371257485, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 5}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '243452', 'license': 'MIT', 'n_tokens_mistral': 1929, 'n_tokens_neox': 1828, 'n_words': 511, 'path': 'Worksheet8/seed_notebook.ipynb', 'repo_name': 'mahesh27dx/molecular_phys', 'size': '5903'}