# Copyright 2012 by the University of Oxford
# Copyright 2014, 2015 by Medical Research Council
# Author: Paul Emsley
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
# 02110-1301, USA

from __future__ import print_function

"""
pyrogen contains tools for converting molecules to mmcif restraints dictionaries,
utilities for retrival, extraction and depction.
"""

import sys
import os
import copy
from subprocess import call
from rdkit import Chem
from rdkit.Chem import AllChem

import coot_git

# Hello from 20170701, library resolution problems?
# $ otool -L .libs/_pyrogen_swig.so -> @rpath substition... how does that work?
import pyrogen_swig as pysw # relies on above rdkit.Chem import AllChem (but ideally should not)
import pyrogen_boost        # ditto
import atom_types

from optparse import OptionParser

import tautomer

import urllib

from jay_util import *

global pyrogen_version
pyrogen_version = "0.0-pre"

global run_mogul
global smiles_dict
run_mogul = False
smiles_dict = False

def make_mogul_ins_file(mogul_ins_file_name, mogul_out_file_name, sdf_file_name):
   f = open(mogul_ins_file_name, 'w')
   if f:
     f.write('mogul molecule file ')
     f.write(sdf_file_name)
     f.write('\n')
     f.write('mogul output   file ')
     f.write(mogul_out_file_name)
     f.write('\n')
     f.write('mogul output distribution all on\n')
     f.write('bond all\n')
     f.write('angle all\n')
#      f.write('torsion all\n')
#      f.write('ring all\n')
     f.write('config output format CSV\n')
     f.write('config output items fragment_type atom_indices query_value nhits mean median sd z-score dmin\n')
     f.write('config search all filter exclude_solvents\n')
     f.write('config output invalid_fragments exclude\n')
     f.close()
   return f


# return True for good, False for bad/not-run
#
def execute_mogul(sdf_file_name, mogul_ins_file_name, mogul_out_file_name):
   f = make_mogul_ins_file(mogul_ins_file_name, mogul_out_file_name, sdf_file_name)
   if f:
      # print 'now run mogul using ins file %s' % mogul_ins_file_name
      if run_mogul:
          state = call(['mogul', '-ins', mogul_ins_file_name])
          return (state == 0)
      else:
          return False
   else:
      return False

def atom_name_from_atomic_number_and_count(element, count, inc):
    name = element
    name += str(count+inc)
    return name

def add_atom_names(mol):
    nz = {}
    atom_names = []
    for atom in mol.GetAtoms():
       try:
          n = atom.GetProp('name')
          atom_names.append(n)
       except KeyError:

          # we want to generate a name that is not already in the atom_names list

          z = atom.GetAtomicNum()
          if z in nz:
             nz[z]  = nz[z] + 1
          else:
             nz[z] = 1;
          ele = atom.GetSymbol().upper()
          # we add inc argument, which gets added to count (nz[z]) in case that we
          # already have a name that matches (previous) return from call to
          # atom_name_from_atomic_number_and_count()
          #
          inc = 0
          name = atom_name_from_atomic_number_and_count(ele, nz[z], inc)
          p_name = pad_atom_name(name, ele)
          # print('c.f.', name, " atom names", atom_names)
          while p_name in atom_names :
             inc += 1
             name = atom_name_from_atomic_number_and_count(ele, nz[z], inc)
             p_name = pad_atom_name(name, ele)
          # print(atom, 'made-name', p_name, ":")

          atom.SetProp("name", p_name)
          atom_names.append(p_name)
    return atom_names

def convert_to_coot_bond_type(rdkit_type):
    out_type = 'single'
    if (rdkit_type == Chem.rdchem.BondType.SINGLE):
        out_type = 'single'
    if (rdkit_type == Chem.rdchem.BondType.AROMATIC):
        out_type = 'aromatic'
    if (rdkit_type == Chem.rdchem.BondType.DOUBLE):
        out_type = 'double'
    if (rdkit_type == Chem.rdchem.BondType.TRIPLE):
        out_type = 'triple'
    if (rdkit_type == Chem.rdchem.BondType.ONEANDAHALF):
        out_type = 'deloc'
    return out_type

def pad_atom_name(name, element):

    padded = name
    if (len(element) == 1):
        if (len(name) == 2):
            padded = ' ' + name + ' '
        if (len(name) == 3):
            padded = ' ' + name
    if (len(element) == 2):
        if (len(name) == 2):
            padded = name + '  '
        if (len(name) == 3):
            padded = name + ' '
    return padded

def is_smiles_file(file_name):
    bits = file_name.rsplit(".")
    if len(bits) > 1:
        return bits[1] == 'smi'
    else:
        return False

def is_comp_id(comp_id):
   return len(comp_id) == 3

def is_mdl_file(file_name):
    bits = file_name.rsplit(".")
    if (len(bits) < 2):
       return False
    else:
       idx = len(bits) - 1
       if (bits[idx] == 'mol'):
          return True
       else:
          if (bits[idx] == 'mdl'):
             return True
          else:
             return False


# return the contents of file_name
def read_file(file_name):
    f = open(file_name)
    r = f.read()
    f.close()
    return r


# return False or a file_name
#
# downloaded file is put in the CCD directory as CCD/x/xyz.cif
#
def get_pdbe_cif_for_comp_id(comp_id):

   try:
      url = 'ftp://ftp.ebi.ac.uk/pub/databases/msd/pdbechem/files/mmcif/' + comp_id + '.cif'
      url = 'https://www.ebi.ac.uk/pdbe/static/files/pdbechem_v2/'        + comp_id + '.cif'
      CCD_dir = 'CCD'
      # file_name = "PDBe-" + comp_id + ".cif"
      first_char = comp_id[0]
      try:
         sub_dir = os.path.join(CCD_dir, first_char)
         file_name = os.path.join(sub_dir, comp_id + ".cif")
         if not os.path.isdir(CCD_dir):
            os.mkdir(CCD_dir)
         if not os.path.isdir(sub_dir):
            os.mkdir(sub_dir)
         if os.path.isfile(file_name):
            return file_name
         else:
            url = 'ftp://ftp.ebi.ac.uk/pub/databases/msd/pdbechem/files/mmcif/' + comp_id + '.cif'
            url = 'https://www.ebi.ac.uk/pdbe/static/files/pdbechem_v2/'        + comp_id + '.cif'
            status = urllib.urlretrieve(url, file_name)
            print('urllib.urllib returned with status', status)
            return file_name

      except OSError as e:
         print(e)
         print("Failed: Can't ftp from", url, "and write file", file_name)

   except IOError as e:
      print(e)
      print("Failed: Can't ftp from", url, "and write file", file_name)

def fetch(comp_id):
    return get_pdbe_cif_for_comp_id(comp_id)

def MolFromFetchedCode(code):
   f = get_pdbe_cif_for_comp_id(code)
   m = pyrogen_boost.MolFromPDBXr(f, code)
   m.Compute2DCoords()
   return m

def make_restraints_for_bond_orders(mol):
    restraints = {}
    bond_list = []
    for bond in mol.GetBonds():
        type = bond.GetBondType()
        coot_bond_type = convert_to_coot_bond_type(type)
        at_1 = bond.GetBeginAtom()
        at_2 = bond.GetEndAtom()
        name_1 = at_1.GetProp('name')
        name_2 = at_2.GetProp('name')
        item = [name_1, name_2, coot_bond_type, 1.0, 1.0]
        bond_list.append(item)
    restraints['_chem_comp_bond'] = bond_list
    restraints['_chem_comp'] = [mol.GetProp('comp_id'),
                                mol.GetProp('comp_id'),
                                mol.GetProp('name'),
                                'non-polymer',
                                mol.GetNumAtoms(),
                                mol.GetNumAtoms(),
                                '.']
    return restraints


# return True if mogul is not run or mogul exe is in place.
# return False if mogul is expected but not found.
def test_for_mogul():
   global run_mogul
   if run_mogul:
      mogol_exe = which('mogul')
      if (mogol_exe == None):
         print("mogul not found in path")
         print("Try: pyrogen --no-mogul")
         return False
      else:
         return True
   else:
      return True # OK, really

# this can throw a TypeError
#
def get_smiles_from_comp_id(comp_id):
   global smiles_dict
   if (not smiles_dict):
      read_smiles_tab('smiles.tab')
   return smiles_dict[comp_id]

# return a dictionary or False (if the file does not exist)
# (can this go inside get_smiles_from_comp_id?)
#
def read_smiles_tab(file_name):
    global smiles_dict
    try:
       smiles_dict = {}
       f = open(file_name)
       lines = f.readlines()
       for line in lines:
           bits = line.rstrip().rsplit()
           smiles_dict[bits[0]] = bits[2]
       f.close()
       return True
    except IOError as e:
       smiles_dict = True # we've tested for it
       return False

# return a pair, the smiles string and the molecule name (which might be blank)
#
def get_smiles_from_file(file_name):
    if not os.path.exists(file_name):
       return False,False
    else:
        f = open(file_name)
        smi_line = f.readline()
        parts = smi_line.split()
        f.close()
        return parts[0], ' '.join(parts[1:])

def make_picture(mol, conf_id, comp_id, output_postfix):

   output_file_name = comp_id + "-" + output_postfix + '.png'
   make_picture_to_file(mol, conf_id, output_file_name)

def make_picture_to_file(mol, conf_id, output_file_name):

   # need to install Pillow (I did it with pip)

   try:
      from rdkit.Chem import Draw
      # The import of Image may change depending on how it was provided.
      # What about pillow? Hmm. Not sure of the details.
      # import Image
      from PIL import Image as Image
      state = Draw.MolToFile(mol, size=(300,300), fileName=output_file_name, confId=conf_id)
      # print 'INFO:: wrote PNG   "' + output_file_name + '"'

      # img = Draw.MolToImage(mol, fitImage=True, size=(900,900))
      # img2 = img.resize((300, 300), Image.ANTIALIAS)
      # img2.save(output_file_name + "resampled.png")

      # testing MolDraw2D code (squiggly mess)
      if True:
         conf_id = AllChem.Compute2DCoords(mol)
         drawer = Draw.MolDraw2DCairo(300,300)
         drawer.DrawMolecule(mol, confId=conf_id)
         drawer.FinishDrawing()
         # c = drawer.GetDrawingText()
   except ImportError as e:
      print('ImportError:', e)
   except ValueError as e:
      print('ValueError in make_picture():', e)

def make_restraints_from_smiles(smiles_string, comp_id, compound_name, mogul_dir, name_stub, pdb_out_file_name, mmcif_dict_name, quartet_planes, quartet_hydrogen_planes, use_mmff, match_atom_names_to_dict_flag, comp_id_list_for_names_match, dict_file_for_names_match):

   if not test_for_mogul():
      # return False
      exit(1)
   m = Chem.MolFromSmiles(smiles_string)
   if compound_name:
       m.SetProp('_Name', compound_name)
   do_hydrogen_atoms_shift = True
   do_minimization = True
   embed_and_optimize = True
   return make_restraints(m, comp_id, mogul_dir, name_stub, pdb_out_file_name, mmcif_dict_name, quartet_planes, quartet_hydrogen_planes, use_mmff, do_minimization, embed_and_optimize, match_atom_names_to_dict_flag, comp_id_list_for_names_match, dict_file_for_names_match, do_hydrogen_atoms_shift)

# return the molecule and return value from make_restraints
#
def make_restraints_from_mdl(mol_file_name, comp_id, mogul_dir, name_stub, pdb_out_file_name, mmcif_dict_name, quartet_planes, quartet_hydrogen_planes, use_mmff, match_atom_names_to_dict_flag, comp_id_list_for_names_match, dict_files_for_names_match):

   if (not (test_for_mogul())):
      # return False, False
      exit(1)

   if not os.path.exists(mol_file_name):
      print("No such file:", mol_file_name)
      exit(1)

   do_hydrogen_atoms_shift = False
   compound_name = '.'
   m = Chem.MolFromMolFile(mol_file_name)
   do_minimization = True
   embed_and_optimize = True
   return m, make_restraints(m, comp_id, mogul_dir, name_stub, pdb_out_file_name, mmcif_dict_name,
                             quartet_planes, quartet_hydrogen_planes, use_mmff, do_minimization, embed_and_optimize,
                             match_atom_names_to_dict_flag, comp_id_list_for_names_match,
                             dict_files_for_names_match, do_hydrogen_atoms_shift)


# return a list of (mol, comp_id) pairs for every ligand in the cif
# file.  Often only one of course.
#
def make_restraints_from_mmcif_dict(cif_file_name_in, comp_id, mogul_dir,
                                    output_dir, output_postfix,
                                    quartet_planes, quartet_hydrogen_planes, use_mmff, do_minimization,
                                    pdb_out_file_name, mmcif_restraints_out_file_name,
                                    do_hydrogen_atoms_shift):

   if not test_for_mogul():
       return [(None, None)]

   if comp_id == "TRY_ALL_COMP_IDS":
      types = pysw.types_from_mmcif_dictionary(cif_file_name_in)
      l = []
      for r_type in types:

         file_name_stub = r_type + "-" + output_postfix
         if options.output_dir != ".":
            file_name_stub = os.path.join(options.output_dir, file_name_stub)

         pdb_out_file_name_local              = file_name_stub + ".pdb"
         mmcif_restraints_out_file_name_local = file_name_stub + ".cif"
         #
         t_mol = make_restraints_from_mmcif_dict_single(cif_file_name_in, r_type, mogul_dir,
                                                        output_postfix,
                                                        quartet_planes,
                                                        quartet_hydrogen_planes, use_mmff, do_minimization,
                                                        pdb_out_file_name_local,
                                                        mmcif_restraints_out_file_name_local,
                                                        do_hydrogen_atoms_shift)
         l.append((t_mol, r_type))
      return l
   else:
       # just the one
       m = make_restraints_from_mmcif_dict_single(cif_file_name_in, comp_id, mogul_dir, output_postfix,
                                                  quartet_planes, quartet_hydrogen_planes, use_mmff, do_minimization,
                                                  pdb_out_file_name, mmcif_restraints_out_file_name,
                                                  do_hydrogen_atoms_shift)
       return [(m, comp_id)]

# return a mol, given a sensible comp_id.
#
# Return None on failure
#
def make_restraints_from_mmcif_dict_single(cif_file_name_in, comp_id, mogul_dir, output_postfix,
                                           quartet_planes, quartet_hydrogen_planes, use_mmff, do_minimization,
                                           pdb_out_file_name, mmcif_restraints_out_file_name,
                                           do_hydrogen_atoms_shift):

   if not test_for_mogul():
       return [(None, None)]

   mogul_file_name_stub = comp_id + '-' + output_postfix # file component of files within mogul_dir

   m = pyrogen_boost.rdkit_mol_chem_comp_pdbx(cif_file_name_in, comp_id)

   if False:  # debugging
      for atom in m.GetAtoms():
         try:
            name    = atom.GetProp('name')
            cipcode = atom.GetProp('_CIPCode')
            print('DEBUG:: make_restraints_from_mmcif_dict_single atom', \
               atom, 'name', name, 'cip-code ', cipcode)
         except KeyError as e:
            print('DEBUG:: pyrogen.py:: make_restraints_from_mmcif_dict_single atom', \
               atom, " with name ", name, ' has no _CIPCode property')
         try:
            r = atom.GetProp('_CIPRank')
            print(atom, "DEBUG:: pyrogen.py:: make_restraints_from_mmcif_dict_single atom cip rank", r)
         except KeyError as e:
            print('DEBUG:: pyrogen.py:: make_restraints_from_mmcif_dict_single atom', \
               atom, "has no _CIPRank")
            pass

   if False:  # debug
      n = m.GetNumConformers()
      print("debug in make_restraints_from_mmcif_dict_single() n conformers:", n)
      if n == 1:
         conf_id = 0
         conf = m.GetConformer(conf_id)
         print(dir(conf))
         n_atoms = conf.GetNumAtoms()
         positions = conf.GetPositions()
         for pos in positions:
            print("debug::   post-rdkit_mol_chem_comp_pdbx()", pos)

   # maybe user didn't select the correct comp_id for the given dictionary mmcif
   if m.GetNumAtoms() == 0:
      print('No atoms for comp_id', comp_id)
      return False
   else :

      name = ''
      try:
         name = m.GetProp('_Name')
      except KeyError:
         print('caught KeyError in make_restraints_from_mmcif_dict_single() trying GetProp _Name')

      embed_and_optimize = False # start from the given coordinates
      # print('############################ make_restraints_from_mmcif_dict_single() calling make_restraints()')
      # print('############################ make_restraints_from_mmcif_dict_single() with args ',
      #       use_mmff, do_minimization, embed_and_optimize)
      return make_restraints(m, comp_id, mogul_dir, mogul_file_name_stub,
                             pdb_out_file_name, mmcif_restraints_out_file_name,
                             quartet_planes, quartet_hydrogen_planes, use_mmff, do_minimization, embed_and_optimize,
                             False, False, False,
                             do_hydrogen_atoms_shift)


def n_hydrogens(mol):
    n_H = 0
    for atom in mol.GetAtoms():
        if atom.GetAtomicNum() == 1:
            n_H += 1
    return n_H

def checked_mkdir(dirname):
    if not os.path.exists(dirname):
       os.makedirs(dirname)
    else:
       if os.path.isdir(dirname):
          pass # this happens most of the time, I imagine
       else:
          print('Stop:: File', dirname, 'exists but is not a directory')


# return sane_H_mol
#
def make_restraints(m, comp_id, mogul_dir, mogul_file_name_stub, pdb_out_file_name, mmcif_dict_name,
                    quartet_planes, quartet_hydrogen_planes, use_mmff, do_minimization, embed_and_optimize,
                    match_atom_names_to_dict_flag,
                    comp_id_list_for_names_match,
                    dict_files_for_names_match, do_hydrogen_atoms_shift):

   # test here (or in calling functions) if m is sane (i.e. is an rdkit molecule)

   if not isinstance(m, Chem.rdchem.Mol):
      print('ERROR:: not a molecule')
      return False

   n_attempts = 20 * m.GetNumAtoms() # default is 10 * number of atoms.

   # pH-dependent protonation or deprotonation
   #
   # do_hydrogen_atoms_shift = True, now passed as an argument, user configuration

   try:
      compound_name = m.GetProp('_Name');
   except KeyError:
      # this happens all the time when we start from a SMILES, users don't need to see it.
      # print 'caught key error in trying to get _Name in make_restraints() for m'
      compound_name = '.'
   except AttributeError as e:
      # Do we need to see this? Perhaps make_restraints() needs to return a status.
      # print 'AttributeError: problem with molecule in make_restraints()', e, ' on object:', m
      return

   m_H = m
   if n_hydrogens(m) == 0:
       m_H = AllChem.AddHs(m)

   if do_hydrogen_atoms_shift:
      # simple sane pH H-exchanges
      sane_H_mol = pyrogen_boost.hydrogen_transformations(m_H)
      # print >>file('sane_H.mol','w+'),Chem.MolToMolBlock(sane_H_mol)
      AllChem.AssignStereochemistry(sane_H_mol, force=True, cleanIt=True)
   else:
      sane_H_mol = m_H

   AllChem.AssignStereochemistry(sane_H_mol);

   n_confs = sane_H_mol.GetNumConformers()
   # print("debug:: there are", n_confs, "conformers in sane_H_mol")

   # conf_id = AllChem.EmbedMolecule(sane_H_mol, maxAttempts=n_attempts)

   if do_minimization:

      if use_mmff:

         if embed_and_optimize:

            conf_id = AllChem.EmbedMolecule(sane_H_mol, AllChem.ETKDG())
            AllChem.MMFFOptimizeMolecule(sane_H_mol, confId=conf_id)

         if False:  # debugging output
            ba = pyrogen_boost.mmff_bonds_and_angles(sane_H_mol) # uses _forcefield_ of the molecule
            n_bonds = ba.bonds_size()
            if n_bonds > 0:
               for i_bond in range(n_bonds):
                  bond = ba.get_bond(i_bond)
                  print(bond.get_idx_1(), bond.get_idx_2(), bond.get_type(), \
                        bond.get_resting_bond_length(), bond.get_sigma())
            n_angles = ba.angles_size()
            if n_angles > 0:
                for i_angle in range(n_angles):
                    angle = ba.get_angle(i_angle)
                    print(angle.get_idx_1(), angle.get_idx_2(), angle.get_idx_3(), \
                          angle.get_resting_angle(), angle.get_sigma())

      else:
         AllChem.UFFOptimizeMolecule(sane_H_mol, confId=conf_id)

   atom_names = add_atom_names(sane_H_mol)
   all_set = atom_types.set_monomer_library_atom_types(sane_H_mol)  # has deloc bonds now, potentially

   # debug sane_H_mol
   if True:
      molblock = Chem.MolToMolBlock(sane_H_mol)
      # print >> file("sane_H_mol.mol",'w'), molblock
      print(molblock, file=open("sane_H_mol.mol",'w'))

   if (all_set != True):
      return False
   else:

      sane_H_mol.SetProp('comp_id', comp_id)
      sane_H_mol.SetProp('name', compound_name)

      sd_local = mogul_file_name_stub + ".sdf"
      sdf_file_name       = os.path.join(mogul_dir, mogul_file_name_stub + '-mogul.sdf')
      mogul_ins_file_name = os.path.join(mogul_dir, mogul_file_name_stub + '-mogul.ins')
      mogul_out_file_name = os.path.join(mogul_dir, mogul_file_name_stub + '-mogul.out')
      Chem.AllChem.ComputeGasteigerCharges(sane_H_mol)

      moguled_mol = pyrogen_boost.mogulify(sane_H_mol) # Nitro bond orders (and other things?)
      if not os.path.isdir(mogul_dir):
          checked_mkdir(mogul_dir)
          if os.path.isdir(mogul_dir):
              mb = Chem.MolToMolBlock(moguled_mol)
              # print >> file(sdf_file_name,'w'), mb
              print(m, file=open(sdf_file_name,'w'))
      else:
          mb = Chem.MolToMolBlock(moguled_mol)
          # print >> file(sdf_file_name,'w'), mb
          print(mb, file=open(sdf_file_name,'w'))


      bor = make_restraints_for_bond_orders(sane_H_mol)

      # print out the set types:
      print_atom_props = False
      if print_atom_props:
          print('--- Atom Props ---')
      for atom in sane_H_mol.GetAtoms():
         charge = atom.GetProp('_GasteigerCharge') # string?
         name   = atom.GetProp('name')
         try:
            atom_type   = atom.GetProp('type_energy')
            is_aromatic = atom.GetIsAromatic()
            hybrid      = atom.GetHybridization()
            f_charge    = float(charge)
            if print_atom_props:
                print("  atom: %s %s type: %s arom: %s hybrid: %s charge: %6.3f" % (name, atom.GetSymbol(),
                                                                                    atom_type.ljust(4),
                                                                                    str(is_aromatic).ljust(5),
                                                                                    str(hybrid).rjust(3),
                                                                                    f_charge))
         except KeyError:
            print("miss", name, atom.GetSymbol(), charge)

      #
      replace_with_mmff_b_a_restraints = False
      if use_mmff:
          replace_with_mmff_b_a_restraints = True

      # execute_mogul() tests if mogul is executable
      #
      mogul_state = execute_mogul(sdf_file_name, mogul_ins_file_name, mogul_out_file_name)
      if mogul_state:

         # Here we need to think about matching to reference
         # dictionary of amino acids (for standard atom names).
         # That function takes a dictionary and a mmdb::Residue.
         # How does that fit in here?
         #
         restraints = pysw.mogul_out_to_mmcif_dict_by_mol(mogul_out_file_name, comp_id,
                                                          compound_name, sane_H_mol, bor,
                                                          mmcif_dict_name, # not used
                                                          quartet_planes,
                                                          quartet_hydrogen_planes,
                                                          replace_with_mmff_b_a_restraints)


         # match_atom_names_to_dict_flag, comp_id_list_for_names_match, dict_file_for_names_match
         if match_atom_names_to_dict_flag:

             restraints = atom_match_dictionary(restraints, sane_H_mol,
                                                comp_id_list_for_names_match,
                                                dict_files_for_names_match)

         pysw.write_restraints(restraints, mmcif_dict_name)
         pysw.regularize_and_write_pdb(sane_H_mol, restraints, comp_id, pdb_out_file_name)

      else:

          # mogul failed or was not in the path:

          if run_mogul == False:

              print('--- "no mogul" path ---')

              # ... but that's OK if we told pyrogen to run without mogul

              # sane_H_mol:
              # print(Chem.MolToMolBlock(sane_H_mol), file=file('debug_sane_H.mol','w+'))

              # debug
              if False:
                 for atom in sane_H_mol.GetAtoms():
                    try:
                       r = atom.GetProp('_CIPRank')
                       print(atom, "DEBUG:: pyrogen.py::make_restraints(): atom CIP rank", r)
                    except KeyError as e:
                       print('DEBUG:: pyrogen.py::make_restraints() atom', atom, "has no _CIPRank")

              do_regularize = True
              if not do_minimization:
                 do_regularize = False

              restraints = pysw.mmcif_dict_from_mol(comp_id, compound_name, sane_H_mol,
                                                    do_minimization,
                                                    mmcif_dict_name,
                                                    quartet_planes, quartet_hydrogen_planes,
                                                    replace_with_mmff_b_a_restraints)

              if restraints == None:
                  print("No restraints")
                  return True # hacked in value

              if match_atom_names_to_dict_flag:

                  restraints = atom_match_dictionary(restraints, sane_H_mol,
                                                     comp_id_list_for_names_match,
                                                     dict_files_for_names_match)
                  pysw.write_restraints(restraints, mmcif_dict_name)

              pysw.write_pdb_from_mol(sane_H_mol, comp_id, pdb_out_file_name)

          else:
              # ... but not if we wanted to use mogul.
              # (We get here if there is a license error for mogul)
              exit(1)

      return sane_H_mol

def atom_match_dictionary(restraints, sane_H_mol, comp_id_list_for_names_match, dict_files_for_names_match):

    template_comp_ids = ['CYS', 'ASP', 'GLU',        'HIS', 'ILE', 'LYS', 'LEU', 'MET',
                         'ASN', 'PRO', 'GLN', 'ARG', 'SER', 'THR', 'VAL', 'TRP', 'TYR',
                         'G',   'A',     'C',   'U',    'GLC', 'MAN']

    if isinstance(comp_id_list_for_names_match, str):
        template_comp_ids = comp_id_list_for_names_match.split(',')

    template_cif_dict_files_names = []
    if isinstance(dict_files_for_names_match, str):
        template_cif_dict_files_names = dict_files_for_names_match.split(',')
        # don't use my set of comp_ids then
        template_comp_ids = []

    success,new_restraints,at_name_list = pysw.match_restraints_to_dictionaries(restraints,
                                                                                template_comp_ids,
                                                                                template_cif_dict_files_names)
    if success:
        for at_name in at_name_list:
            # print at_name
            pass
        n = len(sane_H_mol.GetAtoms())
        if len(restraints['_chem_comp_atom']) == n:
            restraints = new_restraints
            for iat in range(n):
                name = sane_H_mol.GetAtomWithIdx(iat).GetProp('name')
                if name != restraints['_chem_comp_atom'][iat][0]:
                    # print "   changing name from", name, "to", restraints['_chem_comp_atom'][iat][0]
                    sane_H_mol.GetAtomWithIdx(iat).SetProp('name', restraints['_chem_comp_atom'][iat][0]);

    return restraints

def simple_make(mmcif_file_name, comp_id):

   global run_mogul
   run_mogul = False
   mmcif_restraints_out_file_name = comp_id + "-pyrogen.cif"
   pdb_fn = comp_id + "-pyrogen.pdb"
   do_hydrogen_atoms_shift = True
   do_minimization = True
   mol_pairs = make_restraints_from_mmcif_dict(mmcif_file_name,
                                               comp_id,
                                               ".", ".", '.', 'postfix', False, True, True, do_minimization,
                                               pdb_fn, mmcif_restraints_out_file_name, do_hydrogen_atoms_shift)
   for mol_info in mol_pairs:
      (mol, comp_id) = mol_info
   if not mol:
      print('No molecule')


def png_from_mmcif_file(mmcif_file_name_in, comp_id, png_file_name):
   m = pyrogen_boost.rdkit_mol_chem_comp_pdbx(mmcif_file_name_in, comp_id)
   conf_id = 0
   n = m.GetNumConformers()
   if n == 0:
      conf_id = AllChem.Compute2DCoords(m)
   conf = m.GetConformer(conf_id)
   mol_for_drawing = Chem.RemoveHs(m, implicitOnly=False)
   if conf.Is3D():
      # print '3D path'
      conf2D_id = AllChem.Compute2DCoords(mol_for_drawing)
      make_picture_to_file(mol_for_drawing, conf2D_id, png_file_name)
   else:
      # print 'Non-3D path'
      make_picture_to_file(mol_for_drawing, -1, png_file_name)

def coot_png_from_mmcif_file(mmcif_file_name_in, comp_id, png_file_name, n_pixels=300, BackgroundColor=None):
   # change the name of cairo_png_depict
   # BackgroundColor is a hash colour, for example use '#ffffff' for white
   #
   return pyrogen_boost.cairo_png_depict(mmcif_file_name_in, comp_id, png_file_name, n_pixels, BackgroundColor)

def depict(mol, iconf = -1, npx=300, highlightAtoms=[], highlightBonds=None, highlightAtomColours=None, highlightBondColours=None):
    import IPython

    # import Image
    # import ImageFont
    from PIL import Image
    from PIL import ImageFont
    import PIL
    import PIL.ImageDraw
    import io
    try:
       n_confs = mol.GetNumConformers()
       if n_confs == 0:
           iconf = mol.Compute2DCoords()
       s = pyrogen_boost.cairo_png_depict_to_string(mol, iconf, highlightAtoms, highlightBonds, highlightAtomColours, highlightBondColours, npx)
       if len(s) > 0:
          sio = io.BytesIO(s)
          im = Image.open(sio)
          bo = io.BytesIO()
          im.save(bo, 'png')
          IPython.display.display(IPython.display.Image(bo.getvalue()))
       else:
          print('Null image')
    except AttributeError as e:
        # maybe mol was not a RDKit molecule
        print('ERROR::', e)

def coot_depict_to_png_string(mol, iconf=-1, n_px=300, highlightAtoms=[], highlightBonds=None, highlightAtomColours=None, highlightBondColours=None):
       s = pyrogen_boost.cairo_png_depict_to_string(mol, iconf, highlightAtoms, highlightBonds, highlightAtomColours, highlightBondColours, n_px)
       return s


def coot_depict_to_svg_string(mol, iconf=-1, n_px=300, highlightAtoms=[], highlightBonds=None, highlightAtomColours=None, highlightBondColours=None):
       s = pyrogen_boost.cairo_svg_depict_to_string(mol, iconf, highlightAtoms, highlightBonds, highlightAtomColours, highlightBondColours, n_px)
       return s

# make MolFromPDBXr available in pyrogen
def MolFromPDBXr(cif_file_name, comp_id):
    return pyrogen_boost.MolFromPDBXr(cif_file_name, comp_id)

def MolsToGridImage(mols, mols_per_row=3, sub_image_size=(200,200), legends=None):

    # return None on failure to find the font.
    # Locate the ttf font file on disk somewhere and make a truetype
    # font from that (if you can)
    def find_font():

       import distutils
       import distutils.sysconfig
       prfx=distutils.sysconfig.PREFIX

       # this ttf file is hard-coded and copied from dials base/share/fonts
       # we need to install that directory as part of coot and find
       # the file in the installation
       font_file_base_name = "Vera.ttf" # VeraMono
       dir_1 = os.path.join(prfx,  'share')
       dir_2 = os.path.join(dir_1, 'coot')
       dir_3 = os.path.join(dir_2, 'fonts')
       font_file="/usr/share/fonts/truetype/ttf-bitstream-vera/Vera.ttf"
       if os.path.isfile(font_file):
          font=ImageFont.truetype(font_file, 16)
          return font
       else:
          return None

    import IPython
    # import Image
    # import ImageFont
    from PIL import Image
    from PIL import ImageFont
    import PIL
    import PIL.ImageDraw
    import io
    n_rows=(len(mols) // mols_per_row)
    if len(mols) > n_rows*mols_per_row:
       n_rows += 1
    full_size=(sub_image_size[0]*mols_per_row, sub_image_size[1]*n_rows)
    composite=composite = Image.new('RGBA', full_size)
    font=find_font()

    for count,mol in enumerate(mols):
        try:
            # print('mol:', mol)
            s = pyrogen_boost.cairo_png_depict_to_string(mol, -1, None, None, None, None, sub_image_size[0])
            sio = io.BytesIO(s)
            im = Image.open(sio)
            # add label if possible
            try:
               if font:
                  lab = legends[count]
                  text_x_pos = sub_image_size[0]/2 - 10
                  if text_x_pos < 0:
                     text_x_pos = 0
                  text_y_pos = 0.9 * sub_image_size[0]
                  min_to_fit = 15
                  if sub_image_size[1]-text_y_pos < min_to_fit:
                     text_y_pos = sub_image_size[1] - min_to_fit
                  draw = PIL.ImageDraw.Draw(im)
                  # we don't have this attribute (multiline_text). Hmm. Maybe old version of PIL?
                  # draw.multiline_text((text_x_pos, text_y_pos), lab, fill=(10,10,10,255), font=font, align="center")
                  draw.text((text_x_pos, text_y_pos), lab, fill=(10,10,10,255), font=font)
                  # debug print(text_x_pos, text_y_pos)
            except TypeError as e:
               # labels was not a list
               pass
            except NameError as e:
                print(e, "in legends", legends)
            except IOError as e:
               print("something label-related:", e)
            region = im.crop((0,0,sub_image_size[0], sub_image_size[0]))
            i_row = count // mols_per_row
            i_col = count - i_row * mols_per_row
            x_off = i_col * sub_image_size[0]
            y_off = i_row * sub_image_size[0]
            composite.paste(region,(x_off, y_off, x_off+sub_image_size[0], y_off+sub_image_size[0]))
        except ValueError as e:
            print(e, 'count ', count, " for mol", mol)
        except IOError as e:
            print(e, 'count ', count, " for mol", mol)
        except Exception as e:
            print(e, 'count ', count, " for mol", mol)
    b = io.BytesIO()
    composite.save(b, 'png')
    IPython.display.display(IPython.display.Image(b.getvalue()))

def prep(mol):
   m_sh = Chem.RemoveHs(mol)
   m_sh.Compute2DCoords()
   Chem.Kekulize(m_sh)
   return m_sh

def d(mol, npx=300):
   m_sh = prep(mol)
   depict(m_sh, npx=npx)

def score_and_print_tautomers(mol, comp_id, output_postfix, do_drawings):

    results = tautomer.enumerate_tautomers(mol)
    for i in range(len(results)):
        m = results[i]
        s = Chem.MolToSmiles(m)
        print("comp_id :", comp_id, ": SMILES", s, 'score:', tautomer.tautomer_score(m))
        if do_drawings:
            file_name = comp_id + '-tautomer-' + str(i)
            file_name += '-' + options.output_postfix + '.png'
            n = m.GetNumConformers()
            conf_id = 0
            if n == 0:
                conf_id = AllChem.Compute2DCoords(m)
            conf = m.GetConformer(conf_id)
            if conf.Is3D():
                mol_for_drawing = Chem.RemoveHs(m, implicitOnly=False)
                conf2D_id = AllChem.Compute2DCoords(mol_for_drawing)
                make_picture_to_file(mol_for_drawing, conf2D_id, file_name)
            else:
                make_picture_to_file(m, -1, file_name)


if __name__ == "__main__":

    def smiles_and_name_from(smi_raw):
       extension = os.path.splitext(smi_raw)[1]
       smiles_string = ''
       name=''
       if extension == '.smi' or extension == '.smiles':
           if not os.path.exists(smi_raw):
               print("File not found:", smi_raw)
               exit(1)
           else:
               smiles_string,name = get_smiles_from_file(smi_raw)
       else:
         smiles_string = smi_raw
       return smiles_string,name

    parser = OptionParser(usage='pyrogen [options] file-or-SMILES'+
                          '\n       if file-or-SMILES has extension ".smi" or ".smiles" ' +
                          'then it is treated as a file')
    parser.add_option("-c", "--mmcif", dest="mmcif_file_name",
                      help="Make restraints from input mmcif FILE", metavar="FILE")
    parser.add_option("-m", "--mol", dest="sdf_file",
                      help="Make restraints from input sdf/mol FILE", metavar="FILE")
    parser.add_option("-r", "--residue-type", dest="comp_id", default='default',
                      help="Create restraints for this type. Default is LIG")
    parser.add_option("-4", "--quartet-planes", dest="quartet_planes",
                      default=False,
                      help="Use 4-atom plane restraints,\n                    " +
                      "forces --quartet-hydrogens", action="store_true")
    parser.add_option("-H", "--quartet-hydrogens", dest="quartet_hydrogen_planes",
                      default=False,
                      help="Use 4-atom hydrogen plane restraints",
                      action="store_true")
    parser.add_option("-b", "--no-shift-hydrogen-atoms", dest="do_hydrogen_atoms_shift",
                      default=True, help="Stop addition or deletion of Hydrogen atoms for formally-charged atoms",
                      action="store_false")
    parser.add_option("-n", "--no-mogul", dest="use_mogul", default=False, action="store_false",
                      help='Don\'t run CSD Mogul to update bond and angle restraints')
    parser.add_option("-N", '--name', dest='compound_name', default=False,
                      help='Compound name')
    parser.add_option('-S', '--smiles', dest="show_smiles",
                      default=False, action="store_true", help="Write the SMILES for the input molecule")
    parser.add_option("-t", "--tautomers", dest="show_tautomers",
                      default=False, action="store_true",
                      help='Show SMILES for tautomers, don\'t generate restraints')
    parser.add_option("-T", '--tmp-directory', dest='mogul_dir',
                      help='Directory into which the tmp files (e.g. for mogul) are written',
                      default='pyrogen-mogul')
    parser.add_option("-d", '--directory', dest='output_dir',
                      help='Directory into which the output files (e.g. mmCIF and PDB) are written',
                      default='.')
    parser.add_option('-o', '--output-postfix', default='pyrogen',
                      dest='output_postfix',
                      help='string to add to output file names, default is "pyrogen"')
    parser.add_option('-p', '--picture', dest='drawing',
                      help='Additionally output a chemical diagram PNG',
                      action='store_true', default=False)
    parser.add_option('-P', '--preserve-input-coordinates', dest='preserve_input_coordinates',
                      help="Preserve the inputput coordinates (if possible)",
                      action='store_true', default=False)
    parser.add_option('-v', '--version', dest='show_version', default=False,
                      action='store_true', help='Print version information')
    # parser.add_option('-M', '--MMFF', dest='use_mmff', default=False,                       # true now.
    #                  action='store_true', help='Use MMFF fallbacks for bonds and angles')
    parser.add_option('-a', '--no-match-vs-reference-dictionaries', default=False,
                      action='store_true', dest='no_match_names_flag',
                      help="Don't match atom names vs. dictionary molecules (default False)")
    parser.add_option('-R', '--reference-dictionary-files', dest='dict_files_for_names_match',
                      help='Try to match the atom names of the output molecule '+
                      'to this dictionary in these files (comma-separated list)', default=False)
    parser.add_option('-C', '--reference-dictionary-comp-ids', dest='comp_id_list_for_names_match',
                      help='Try to match the atom names of the output molecule to these comp-ids' +
                      ' (comma-separated list)',
                      default=False)
    parser.add_option('-w', '--wwPDB', default=False, dest="wwPDB", action="store_true",
                      help='Fetch the wwPDB ligand definition and use that')
    parser.add_option('-f', dest="fetch", default=False,
                      help='(Just) fetch from the PDBe the CCD entry for the given compound-id')
    parser.add_option("-q", "--quiet",
                      action="store_false", dest="verbose", default=True,
                      help="print less messages")

    (options, args) = parser.parse_args()
    # print('DEBUG:: options:', options)
    # print('DEBUG:: args:', args)

    if len(sys.argv) == 0:
       print("Usage: pyrogen --help")

    if options.show_version:
       print('pyrogen-' + pyrogen_version, "revision", coot_git.revision_count())

    if options.fetch:
        fetch(options.fetch)

    comp_id = options.comp_id
    if options.comp_id == 'default':
        comp_id = 'LIG'
    if options.mmcif_file_name != None:
        if options.comp_id == 'default':
            comp_id = 'TRY_ALL_COMP_IDS'

    file_name_stub = comp_id + '-' + options.output_postfix

    if options.output_dir != ".":
       file_name_stub = os.path.join(options.output_dir, file_name_stub)

    pdb_out_file_name              = file_name_stub + '.pdb'
    mmcif_restraints_out_file_name = file_name_stub + '.cif'

    # this is a bit ugly, perhaps. this value is inspected inside
    # the following functions
    # (run_mogul is a global, default no)
    if options.use_mogul == True:
       run_mogul = True
    if run_mogul:
       if len(options.mogul_dir) > 0:
          if options.mogul_dir[0] == '-':
             print('Stop:: you probably didn\'t mean that you wanted',options.mogul_dir, 'as your tmp directory.')
             exit(1)
          checked_mkdir(options.mogul_dir)

    if options.show_tautomers or options.show_smiles:

        # ------------------------ Tautomers and SMILES  ---------------------------------------------

        mol = False
        if len(args) > 0:
            smi_raw = args[0]
            smiles,compound_name = smiles_and_name_from(smi_raw)
            mol = Chem.MolFromSmiles(smiles)
        else:
            if options.sdf_file != None:
                mol = Chem.MolFromMolFile(options.sdf_file)
            else:
                if options.mmcif_file_name != None:
                    types = pysw.types_from_mmcif_dictionary(options.mmcif_file_name)
                    print('-- tautomer mode: mmcif file types:', types)
                    for type in types:
                        mol_local = pyrogen_boost.rdkit_mol_chem_comp_pdbx(options.mmcif_file_name, type)
                        score_and_print_tautomers(mol_local, type, options.output_postfix, options.drawing)

        if mol:
           if options.show_tautomers:
              score_and_print_tautomers(mol, comp_id, options.output_postfix, options.drawing)
           if options.show_smiles:
              s = Chem.MolToSmiles(mol);
              print(s)

    else:

        # ------------------------ dict-build-mode ---------------------------------------------------

        mmcif_file_name = options.mmcif_file_name
        use_mmff = True
        do_minimization = True
        if options.preserve_input_coordinates == True:
           do_minimization = False
        # shall we go get the dictionary?
        if options.wwPDB:
           mmcif_file_name = get_pdbe_cif_for_comp_id(comp_id)
           if os.path.isfile(mmcif_file_name):
              pass # good
           else:
              print("Missing downloaded file for comp-id:",  comp_id)
              exit(2)

        # JED mode for hydrogen planes
        #
        quartet_hydrogen_planes = options.quartet_hydrogen_planes
        if options.quartet_planes:
            quartet_hydrogen_planes = True

        match_names_flag = True
        if options.no_match_names_flag:
            match_names_flag = False


        if mmcif_file_name:
            mol_pairs = make_restraints_from_mmcif_dict(mmcif_file_name,
                                                        comp_id,
                                                        options.mogul_dir,
                                                        options.output_dir,
                                                        options.output_postfix,
                                                        options.quartet_planes,
                                                        quartet_hydrogen_planes,
                                                        use_mmff, do_minimization,
                                                        pdb_out_file_name,
                                                        mmcif_restraints_out_file_name,
                                                        options.do_hydrogen_atoms_shift)

            # this needs to be in a try block, I suppose, for example if the mmcif file
            # does not exist.

            for mol_info in mol_pairs:
                (mol, comp_id) = mol_info
                if not mol:
                    print('No molecule')
                else:
                    # Happy path
                    if options.drawing:
                        # make_picture() by default draws the first conformer in the given molecule.
                        # For mol, that is a 3D conformer.  We want to draw a nice 2D diagram
                        #
                        mol_for_drawing = Chem.RemoveHs(mol, implicitOnly=False)
                        conf2D_id = AllChem.Compute2DCoords(mol_for_drawing)
                        conf = mol.GetConformer(conf2D_id)
                        Chem.WedgeMolBonds(mol_for_drawing, conf)
                        make_picture(mol_for_drawing, conf2D_id, comp_id, options.output_postfix)

        else:

            if options.sdf_file != None:
                (mol, results) = make_restraints_from_mdl(options.sdf_file, comp_id,
                                                          options.mogul_dir, file_name_stub,
                                                          pdb_out_file_name,
                                                          mmcif_restraints_out_file_name,
                                                          options.quartet_planes,
                                                          quartet_hydrogen_planes,
                                                          use_mmff,
                                                          match_names_flag,
                                                          options.comp_id_list_for_names_match,
                                                          options.dict_files_for_names_match)
                if options.drawing:
                    make_picture(mol, -1, comp_id, options.output_postfix)

            else:

                if len(args) > 0:
                    smi_raw = args[0]
                    smiles,compound_name_from_file = smiles_and_name_from(smi_raw)
                    compound_name=False
                    if len(compound_name_from_file) > 0:
                       compound_name = compound_name_from_file
                    print("#####", compound_name, options.compound_name)
                    if isinstance(options.compound_name, str):
                         compound_name = options.compound_name
                    status = make_restraints_from_smiles(smiles, comp_id, compound_name,
                                                         options.mogul_dir, file_name_stub,
                                                         pdb_out_file_name,
                                                         mmcif_restraints_out_file_name,
                                                         options.quartet_planes,
                                                         quartet_hydrogen_planes,
                                                         use_mmff,
                                                         match_names_flag,
                                                         options.comp_id_list_for_names_match,
                                                         options.dict_files_for_names_match)
                    if options.drawing:
                        mol = Chem.MolFromSmiles(smiles)
                        make_picture(mol, -1, comp_id, options.output_postfix)