Mutate model

The script below takes a given PDB file, and mutates a single residue. The residue's position is then optimized, and the unoptimized and optimized energies are reported.

#!python
import sys
import os

#
#  mutate_model.py
#    
#     Usage:   modCVS - < mutate_model.py modelname respos resname chain > logfile
#    
#     Example: modCVS - < mutate_model.py 1t29 1699 LEU A > 1t29.log
#
#
#  Creates a single in silico point mutation to sidechain type and at residue position
#  input by the user, in the structure whose file is modelname.pdb
#  The conformation of the mutant sidechain is optimized by conjugate gradient and
#  refined using some MD.
#


def optimize(mdl1):
    #conjugate gradient
    for step in range(0,len(mdl1.schedule)):
        mdl1.schedule.step = step + 1
        mdl1.optimize(optimization_method=1,max_iterations=200,min_atom_shift=0.001)
    #md
    refine(mdl1)
    mdl1.optimize(optimization_method=1,max_iterations=200,min_atom_shift=0.001)



#molecular dynamics
def refine(mdl1):
    # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
    init_vel = True
    for (its, equil, temps) in ((200, 20, (150.0, 250.0, 400.0, 700.0, 1000.0)),
                                (200, 600, (1000.0, 800.0, 600.0, 500.0, 400.0, 300.0))):
        for temp in temps:
            mdl1.optimize(cap_atom_shift=0.39, md_time_step=4.0, optimization_method=3,
                          md_return='FINAL', init_velocities=init_vel, temperature=temp,
                          max_iterations=its, equilibrate=equil)
            init_vel = False



#use homologs and dihedral library for dihedral angle restraints
def make_restraints(mdl1, aln):
   rsr = mdl1.restraints
   rsr.clear()
   for typ in ('stereo', 'phi-psi_binormal'):
       rsr.make(restraint_type=typ, aln=aln, spline_on_site=True)
   for typ in ('omega', 'chi1', 'chi2', 'chi3', 'chi4'):
       rsr.make(restraint_type=typ+'_dihedral', spline_range=4.0, spline_dx=0.3,
                spline_min_points = 5, aln=aln, spline_on_site=True)

#first argument
modelname, respos, restyp, chain, = sys.argv[1:]


log.verbose()

env = environ(rand_seed=-49837)
env.io.hetatm = True
#soft sphere potential
env.edat.dynamic_sphere=False
#lennard-jones potential (more accuate)
env.edat.dynamic_lennard=True
env.edat.contact_shell = 4.0
env.edat.update_dynamic = 0.39

# Read customized topology file with phosphoserines (or standard one)
env.libs.topology.read(file='$(LIB)/top_heav.lib')

# Read customized CHARMM parameter library with phosphoserines (or standard one)
env.libs.parameters.read(file='$(LIB)/par.lib')


# Read the original PDB file and copy its sequence to the alignment array:
mdl1 = model(env, file=modelname)
ali = alignment(env)
ali.append_model(mdl1, atom_files=modelname, align_codes=modelname)

#set up the mutate residue selection segment
mdl1.pick_atoms(pick_atoms_set=1, selection_from ='ALL', res_types ='ALL',
                selection_segment = (respos+':'+chain, respos+':'+chain),
                atom_types = 'ALL', selection_status = 'INITIALIZE')

#perform the mutate residue operation
mdl1.mutate(residue_type=restyp)

#get two copies of the sequence.  A modeller trick to get things set up
ali.append_model(mdl1, align_codes=modelname)

# Generate molecular topology for mutant
mdl1.generate_topology(ali, sequence=modelname)


# Transfer all the coordinates you can from the template native structure
# to the mutant (this works even if the order of atoms in the native PDB 
# file is not standard):
#here we are generating the model by reading the template coordinates
mdl1.transfer_xyz(ali)

# Build the remaining unknown coordinates
mdl1.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')

#yes model2 is the same file as model1.  It's a modeller trick.
mdl2 = model(env, file=modelname)

#required to do a transfer_res_numb
#ali.append_model(mdl2, atom_files=modelname, align_codes=modelname)
#transfers from "model 2" to "model 1"
mdl1.res_num_from(mdl2,ali)

#It is usually necessary to write the mutated sequence out and read it in
#before proceeding, because not all sequence related information about MODEL
#is changed by this command (e.g., internal coordinates, charges, and atom
#types and radii are not updated).

mdl1.write(file=modelname+restyp+respos+'.tmp')
mdl1.read(file=modelname+restyp+respos+'.tmp')

#set up restraints before computing energy
#we do this a second time because the model has been written out and read in,
#clearing the previously set restraints
make_restraints(mdl1, ali)

#a non-bonded pair has to have at least as many selected atoms
mdl1.env.edat.nonbonded_sel_atoms=1

mdl1.schedule.make(library_schedule=6)

#only optimize the selected residue (in first pass, just atoms in selected
#residue, in second pass, include nonbonded neighboring atoms)
#set up the mutate residue selection segment
mdl1.pick_atoms(pick_atoms_set=1, selection_from ='ALL', res_types ='ALL',
                selection_segment = (respos+':'+chain, respos+':'+chain),
                atom_types = 'ALL', selection_status = 'INITIALIZE')

mdl1.restraints.unpick_all()
mdl1.restraints.pick()

mdl1.energy()

mdl1.randomize_xyz(deviation=4.0)

mdl1.env.edat.nonbonded_sel_atoms=2
optimize(mdl1)

#feels environment (energy computed on pairs that have at least one member
#in the selected)
mdl1.env.edat.nonbonded_sel_atoms=1
optimize(mdl1)

mdl1.energy()

#give a proper name
mdl1.write(file=modelname+restyp+respos+'.pdb')

#delete the temporary file
os.remove(modelname+restyp+respos+'.tmp');