0

I am running a model evaluation protocol for Modeller. It evaluates every model and writes its result to a separate file. However I have to run it for every model and write to a single file.

This is the original code:

from modeller import *
from modeller.scripts import complete_pdb

log.verbose()    # request verbose output
env = environ()
env.libs.topology.read(file='$(LIB)/top_heav.lib') # read topology
env.libs.parameters.read(file='$(LIB)/par.lib') # read parameters

# read model file
mdl = complete_pdb(env, 'TvLDH.B99990001.pdb')

# Assess all atoms with DOPE:
s = selection(mdl)
s.assess_dope(output='ENERGY_PROFILE NO_REPORT', file='TvLDH.profile',
              normalize_profile=True, smoothing_window=15)

I added a loop to evaluate every model in a single run, however I am creating several files (one for each model) and I want is to print all evaluations in a single file

from modeller import *
from modeller.scripts import complete_pdb

log.verbose()    # request verbose output
env = environ()
env.libs.topology.read(file='$(LIB)/top_heav.lib') # read topology
env.libs.parameters.read(file='$(LIB)/par.lib') # read parameters

#My loop starts here
for i in range (1,1001):
    number=str(i)
    if i<10:
        name='000'+number
    else:
        if i<100:
           name='00'+number
        else:
             if i<1000:
                name='0'+number
             else:
                  name='1000'

# read model file
mdl = complete_pdb(env, 'TcP5CDH.B9999'+name+'.pdb')

# Assess all atoms with DOPE: this is the assesment that i want to print in the same file
s = selection(mdl)
savename='TcP5CDH.B9999'+name+'.profile'
s.assess_dope(output='ENERGY_PROFILE NO_REPORT', 
              file=savename,
              normalize_profile=True, smoothing_window=15)

As I am new to programming, any help will be very helpful!

JoeG
  • 7,191
  • 10
  • 60
  • 105

1 Answers1

0

Welcome :-) Looks like you're very close. Let's introduce you to using a python function and the .format() statement.

Your original has a comment line # read model file, which looks like it could be a function, so let's try that. It could look something like this.

from modeller import *
from modeller.scripts import complete_pdb

log.verbose()    # request verbose output
# I'm assuming this can be done just once
# and re-used for all your model files...
# (if not, the env stuff should go inside the
# read_model_file() function.
env = environ()
env.libs.topology.read(file='$(LIB)/top_heav.lib') # read topology
env.libs.parameters.read(file='$(LIB)/par.lib') # read parameters

def read_model_file(file_name):
    print('--- read_model_file(file_name='+file_name+') ---')
    mdl = complete_pdb(env, file_name)
    # Assess all atoms with DOPE:
    s = selection(mdl)
    output_file = file_name+'.profile'
    s.assess_dope(
        output='ENERGY_PROFILE NO_REPORT',
        file=output_file,
        normalize_profile=True,
        smoothing_window=15)

for i in range(1,1001):
    file_name = 'TcP5CDH.B9999{:04d}pdb'.format(i)
    read_model_file(file_name)

Using .format() we can get ride of the multiple if-statement checks for 10? 100? 1000? Basically .format() replaces {} curly braces with the argument(s). It can be pretty complex but you don't need to digetst all of it.

Example: 'Hello {}!'.format('world') yields Hello world!. The {:04d} stuff uses formatting, basically that says "Please make a 4-character wide digit-substring and zero-fill it, so you should get '0001', ..., '0999', '1000'. Just {:4d} (no leading zero) would give you space padded results (e.g. ' 1', ..., ' 999', '1000'. Here's a little more on the zero-fill: Display number with leading zeros

jgreve
  • 1,225
  • 12
  • 17