from __future__ import division
import iotbx.pdb
from scitbx.array_family import flex

pdb_str="""\
CRYST1   17.781   14.493   17.475  90.00  90.00  90.00 P 1
ATOM      1  N   GLY A   1       5.043   9.493  11.271  1.00 16.77           N
ATOM      2  CA  GLY A   1       5.000   9.088   9.820  1.00 16.57           C
ATOM      3  C   GLY A   1       6.037   8.021   9.588  1.00 16.16           C
ATOM      4  O   GLY A   1       6.529   7.402  10.550  1.00 16.78           O
ATOM      5  N   ASN A   2       6.396   7.804   8.324  1.00 15.02           N
ATOM      6  CA  ASN A   2       7.530   6.919   8.000  1.00 14.10           C
ATOM      7  C   ASN A   2       8.811   7.418   8.596  1.00 13.13           C
ATOM      8  O   ASN A   2       9.074   8.623   8.595  1.00 11.91           O
ATOM      9  CB  ASN A   2       7.706   6.762   6.510  1.00 15.38           C
ATOM     10  CG  ASN A   2       6.468   6.223   5.861  1.00 14.08           C
ATOM     11  OD1 ASN A   2       6.027   5.108   6.185  1.00 17.46           O
ATOM     12  ND2 ASN A   2       5.848   7.036   5.000  1.00 11.72           N
ATOM     13  N   ASN A   3       9.614   6.471   9.074  1.00 12.26           N
ATOM     14  CA  ASN A   3      10.859   6.785   9.758  1.00 11.74           C
ATOM     15  C   ASN A   3      12.097   6.213   9.064  1.00 11.10           C
ATOM     16  O   ASN A   3      12.180   5.000   8.817  1.00 10.42           O
ATOM     17  CB  ASN A   3      10.793   6.259  11.211  1.00 12.15           C
ATOM     18  CG  ASN A   3      12.046   6.620  12.030  1.00 12.82           C
ATOM     19  OD1 ASN A   3      12.350   7.806  12.241  1.00 15.05           O
ATOM     20  ND2 ASN A   3      12.781   5.596  12.475  1.00 13.48           N
"""

def r_factor(x,y):
  x = abs(x).data()
  y = abs(y).data()
  return flex.sum(flex.abs(x-y))/flex.sum(flex.abs(x+y))*200.
  
def run():
  pdb_inp = iotbx.pdb.input(source_info=None, lines=pdb_str)
  xray_structure = pdb_inp.xray_structure_simple()
  for d_min in [0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0]:
    xray_structure.scattering_type_registry(d_min=d_min, table="it1992")
    fc_direct = xray_structure.structure_factors(d_min=d_min, 
      algorithm="direct").f_calc()
    fc_fft = fc_direct.structure_factors_from_scatterers(
      xray_structure=xray_structure, 
      algorithm="fft").f_calc()
    print "d_min: %3.1f r_factor=%6.4f"%(d_min, r_factor(fc_direct, fc_fft))
   
if __name__=='__main__':
  run()