build_dimer public version (includes BioPython, PeptideBuilder dependencies)

build_dimer.py

# %%file build_dimer.py
import os
from PeptideBuilder import Geometry
import PeptideBuilder
import Bio.PDB
from simtk.unit import *
from simtk.openmm.app import *
import numpy as np

def buildPeptide(peptide, customAngles=False, phi=-40, psi=60, omega=180):
    """
    construct a peptide sequence pdb file
    :param peptide:
    :return None:
    """
    geo = Geometry.geometry(peptide[0])
    
    if customAngles:
        geo.phi = phi
        geo.psi = psi
        geo.omega = omega
    
    structure = PeptideBuilder.initialize_res(peptide[0])
    
    for i in range(1, len(peptide)):
        geo = Geometry.geometry(peptide[i])
        
        if customAngles:
            geo.phi = phi
            geo.psi = psi
            geo.omega = omega
        
        PeptideBuilder.add_residue(structure, geo)
            
    PeptideBuilder.add_terminal_OXT(structure) # OpenMM will not run without this, but LightDock will not run with it. Solution, add terminal oxygen in prepPDB after docking

    out = Bio.PDB.PDBIO()
    out.set_structure(structure)
    
    if customAngles:
        filename = 'peptide_{}_customAngles.pdb'.format(peptide)
    else:
        filename = 'peptide_{}.pdb'.format(peptide)
    
    out.save(filename)
    
    return filename

def addH(structure, pH):
    """
    protonate a given structure @ a certain pH
    :param file:
    :return:
    """
    pdb = PDBFile(structure)
    topology = pdb.topology
    positions = pdb.positions
    modeller = Modeller(topology, positions)
    modeller.addHydrogens(pH=pH)
    PDBFile.writeFile(modeller.topology, modeller.positions, open(structure.split('.')[0] + '_H.pdb', 'w'))
    
    return structure.split('.')[0] + '_H.pdb'

def runAll(peptide, customAngles=False, minimization=True):
    filename = buildPeptide(peptide, customAngles=customAngles)
    file_to_simul = addH(filename, 7.4) # will create a file called peptide_H.pdb -> this is the one we shall use 

    ## Simulating the peptide dimer - final pdb file will be called "finished_run.pdb"
    pdb = PDBFile(file_to_simul)
    forcefield = ForceField('amber14-all.xml', 'tip3p.xml')
    system = forcefield.createSystem(pdb.topology, nonbondedMethod=NoCutoff)

    K = int(1e4) # K = 10,000 needed for good constraints, K = 100,000 blows up

    force = CustomTorsionForce('0.5*K*dtheta^2; dtheta = min(diff, 2*' 
                               + str(round(np.pi, 3)) + '-diff); diff = abs(theta - theta0)')
    
    force.addGlobalParameter('K', K)
    force.addPerTorsionParameter('theta0')

    force.addTorsion(0, 4, 7, 9, (60 * np.pi / 180,))
    force.addTorsion(4, 7, 9, 11, (180 * np.pi / 180,))
    system.addForce(force)

    integrator = LangevinIntegrator(310*kelvin, 1/picosecond, 0.002*picosecond)
    simulation = Simulation(pdb.topology, system, integrator)
    simulation.context.setPositions(pdb.positions)

    if minimization:
        simulation.minimizeEnergy()

    if customAngles:
        if minimization:
            simulation.reporters.append(PDBReporter('finished_run_wForce{}_customAngles_minim_peptide_{}.pdb'.format(K, peptide), 100))
        else:
            simulation.reporters.append(PDBReporter('finished_run_wForce{}_customAngles_peptide_{}.pdb'.format(K, peptide), 100))
    else:
        if minimization:
            simulation.reporters.append(PDBReporter('finished_run_wForce{}_minim_peptide_{}.pdb'.format(K, peptide), 100))
        else:
            simulation.reporters.append(PDBReporter('finished_run_wForce{}_peptide_{}.pdb'.format(K, peptide), 100))

    simulation.step(5000)

    return

peptide = 'GG'

# for i in tqdm(range(100)):
runAll(peptide, customAngles=True)