from __future__ import print_function
from modeller import *
from modeller.scripts import complete_pdb

log.verbose()
env = Environ()
env.libs.topology.read('${LIB}/top_heav.lib')
env.libs.parameters.read('${LIB}/par.lib')

map='../step6_fit/TvLDH.10A.mrc'  
resolution=10.0  
box_size=48 
apix=1.88
x=-26.742; y=-9.5205; z=-10.375 #origin

m1 = Model(env, file='../step6_fit/TvLD_1_1.pdb')
sel = Selection(m1).only_atom_types('CA')
scal = physical.Values(default=0., em_density=1.0)
results = []
filenames = ['../step6_fit/TvLD_1_1.pdb']
for i in range(1, 6):
    filenames.append('TvLDH.BL%04d0001.pdb' % i)
for fname in filenames:
    m2 = complete_pdb(env, fname)
    aln = Alignment(env)
    aln.append_model(m1, align_codes='fitted')
    aln.append_model(m2, align_codes='model')
    sel.superpose(m2, aln)
    den = Density(env, file=map, em_density_format='MRC',
                  voxel_size=apix, resolution=resolution, em_map_size=box_size,
                  density_type='GAUSS', px=x,py=y,pz=z)
    m2.env.edat.density = den
    molpdf, terms = Selection(m2).energy(schedule_scale=scal)
    results.append((fname, -molpdf))

print("\n%-40s %12s" % ('Filename', 'CCF'))
print("\n".join(['%-40s %12.4f' % x for x in results]))
