#  ============================================================================
#
#  This script calculates the energy of a specific "fit" of an atomic structure
#  in a cryo-EM density map using MODELLER/Mod-EM (CCF, stereo-chemical
#  and non-bonded interactions terms).
#
#  =======================  Maya Topf, 19 Dec 2008 =============================

from modeller import *
import shutil
import sys, os, os.path
import string

class energy:

   def __init__ (self, path, code, run_num, rand, em_map_file, input_pdb_file, format, apix, box_size, res, x, y, z):
      env = environ(rand_seed=-rand) 
      env.libs.topology.read(file='$(LIB)/top_heav.lib')
      env.libs.parameters.read(file='$(LIB)/par.lib')
      log.verbose()

      den = density(env, file=em_map_file, em_density_format=format,
                   em_map_size=box_size, density_type='GAUSS', voxel_size=apix,
                   resolution=res,px=x,py=y,pz=z)   
      
      env.edat.density = den
      env.edat.dynamic_sphere = True
      
# read pdb
      aln = alignment(env)
      mdl2 = model(env, file=input_pdb_file)

      aln.append_model(mdl2, align_codes=input_pdb_file,
                      atom_files=input_pdb_file)
      mdl = model(env, file=input_pdb_file)
      aln.append_model(mdl, align_codes=code, atom_files=code)
      aln.align(gap_penalties_1d=(-600, -400))
      mdl.clear_topology()
      mdl.generate_topology(aln[input_pdb_file])
      mdl.transfer_xyz(aln)
      mdl.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')

# Build restraints
      sel_all = selection(mdl)
      mdl.restraints.make(sel_all, restraint_type='stereo', spline_on_site=False)

      mdl.restraints.make(sel_all, aln=aln, restraint_type='phi-psi_binormal',
                 spline_on_site=True,
                 residue_span_range=(0, 99999))

# evaluate energy of the initial model (last MD model)
      scal = physical.values(default=1.0, em_density=10000)   
      print "energy all (before CG minimization)" 
      print "prot name: %s" %input_pdb_file
      eval = sel_all.energy(schedule_scale=scal)