coots.py

import os
import math
import subprocess
from dimple import utils

M_SQRT1_2 = 0.5**0.5

def find_path():
    if os.name == 'nt':
        for path in ['C:/WinCoot/wincoot.bat',     # since WinCoot 0.8.8
                     'C:/WinCoot/runwincoot.bat',  # WinCoot prior to 0.8.8
                     utils.cbin('coot.bat')]:      # CCP4 script added in 2018
            if os.path.exists(path):
                return path
        utils.comment('\nNote: WinCoot not found.\n')
    else:
        return utils.syspath('coot')

# returns version string
def find_version(coot_path):
    if not coot_path:
        return None
    # On Windows reading output from runwincoot.bat is not reliable.
    # "C:/Windows/bin/coot-real.exe --version" used to work, but in 8.1
    # exe was moved to libexec/coot-bin.exe and DLLs are in different dir
    if os.name == 'nt':
        return 'WinCoot, with python'
    try:
        return subprocess.check_output([coot_path, '--version']).decode()
    except (subprocess.CalledProcessError, OSError):
        return None

def basic_script(pdb, refl, normal_map, center, toward, white_bg):
    same_dir = (os.path.dirname(pdb) == '' and os.path.dirname(refl) == '')
    lines = ['#!/usr/bin/env coot',
             '# python script for coot - generated by dimple',
             'set_nomenclature_errors_on_read("ignore")']
    if same_dir:
        lines += ['import inspect, os',
                  'this_file = inspect.getfile(inspect.currentframe())',
                  'this_dir = os.path.dirname(os.path.abspath(this_file))']
    if white_bg:
        lines += ['set_background_colour(1.0, 1.0, 1.0)']
    if pdb:
        if same_dir:
            lines += ['molecule = read_pdb(os.path.join(this_dir, "%s"))' % pdb]
        else:
            lines += ['molecule = read_pdb("%s")' % pdb]
    if center:
        lines += ['set_rotation_centre(%g, %g, %g)' % center,
                  'set_zoom(30.)']
        if toward:
            lines += ['set_view_quaternion(%g, %g, %g, %g)'
                      % view_as_quat(center, toward)]
    if same_dir:
        lines += ['refl = os.path.join(this_dir, "%s")' % refl]
    else:
        lines += ['refl = "%s"' % refl]

    if normal_map:
        lines += [
            'map21 = make_and_draw_map(refl, "FWT", "PHWT", "", 0, 0)',
            'map11 = make_and_draw_map(refl, "DELFWT", "PHDELWT", "", 0, 1)']
    else:
        lines += [
            'm = read_phs_and_make_map_using_cell_symm_from_previous_mol(refl)',
            'set_last_map_contour_level_by_sigma(4)']

    return '\n'.join(lines) + '\n'

# view from p1 to p2, to be passed to set_view_quaternion()
def view_as_quat(p1, p2):
    if p1 is None or p2 is None:
        return (0., 0., 0., 1.)
    d = (p2[0]-p1[0], p2[1]-p1[1], p2[2]-p1[2])
    length = math.sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2])
    d = (d[0]/length, d[1]/length, d[2]/length)
    # ref vector: 0 0 -1 (?)
    # cross product (a2 b3 - a3 b2, ..., ...)
    prod = (d[1], -d[0], 0)
    quat = (prod[0], prod[1], prod[2], 1-d[2])
    qlen = math.sqrt(sum(a*a for a in quat))
    return (quat[0]/qlen, quat[1]/qlen, quat[2]/qlen, quat[3]/qlen)


def mult_quat(q1, q2):
    x, y, z, w = q1
    ax, ay, az, aw = q2
    return (w*ax + x*aw + y*az - z*ay,
            w*ay + y*aw + z*ax - x*az,
            w*az + z*aw + x*ay - y*ax,
            w*aw - x*ax - y*ay - z*az)


def r3d_script(center, toward, blobname):
    #quat0 = (0., 0., 0., 1.)
    quat0 = view_as_quat(center, toward)
    quaternions = [quat0,
                   mult_quat(quat0, (0., M_SQRT1_2, 0., M_SQRT1_2)),
                   mult_quat(quat0, (M_SQRT1_2, 0., 0., M_SQRT1_2))]

    # Coot function raster3d() creates file.r3d, make_image_raster3d() also
    # calls render program and opens image (convenient for testing)
    basenames = []
    script = """
set_rotation_centre(%g, %g, %g)
set_zoom(30.)""" % center
    for n, quat in enumerate(quaternions):
        script += """
set_view_quaternion(%g, %g, %g, %g)""" % quat
        basename = '%sv%d' % (blobname, n+1)
        script += """
graphics_draw() # this is needed only for coot in --no-graphics mode
raster3d("%s.r3d")""" % basename
        basenames.append(basename)
    return script, basenames