from __future__ import division
import iotbx.pdb
import mmtbx.model
import mmtbx.f_model
from iotbx import reflection_file_reader
from scitbx.array_family import flex

def run():
  pdb_inp = iotbx.pdb.input(file_name = "1yjp.pdb")
  model = mmtbx.model.manager(model_input = pdb_inp)
  miller_arrays = reflection_file_reader.any_reflection_file(file_name =
    "1yjp.mtz").as_miller_arrays()
  for ma in miller_arrays:
    print(ma.info().label_string())
    if(ma.info().label_string()=="FOBS,SIGFOBS"):
      f_obs = ma
  fmodel = mmtbx.f_model.manager(
    f_obs          = f_obs,
    xray_structure = model.get_xray_structure(),
    target_name    = 'ls_wunit_k1')
  fmodel.update_all_scales()
  print("r_work=%6.4f r_free=%6.4f"%(fmodel.r_work(), fmodel.r_free()))
  tf = fmodel.target_functor()(compute_gradients=True)
  print(list(tf.d_target_d_site_cart()))
  # Your way:
  print("-"*79)
  xrs = fmodel.xray_structure
  xrs.scatterers().flags_set_grads(state=False)
  xrs.scatterers().flags_set_grad_site(
    iselection = xrs.all_selection().iselection())
  # This is same as above
  #from cctbx import xray
  #xray.set_scatterer_grad_flags(
  #    scatterers=xrs.scatterers(),
  #    site=True)
  gradients = flex.vec3_double(
    fmodel.one_time_gradients_wrt_atomic_parameters(site=True).packed())
  print(list(gradients))
  

if(__name__ == "__main__"):
  run()