File: crd.cpp

package info (click to toggle)
gemmi 0.6.5%2Bds-3
  • links: PTS, VCS
  • area: main
  • in suites: trixie
  • size: 5,836 kB
  • sloc: cpp: 54,719; python: 4,743; ansic: 3,972; sh: 384; makefile: 73; f90: 42; javascript: 12
file content (150 lines) | stat: -rw-r--r-- 5,609 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
// Copyright 2017-2022 Global Phasing Ltd.

#include <stdio.h>             // for fprintf, stderr, putc
#include <iostream>            // for cerr
#include <exception>           // for exception
#include "gemmi/crd.hpp"       // for prepare_refmac_crd
#include "gemmi/to_cif.hpp"    // for write_cif_block_to_stream
#include "gemmi/fstream.hpp"   // for Ofstream
#include "gemmi/monlib.hpp"    // for MonLib, read_monomer_lib
#include "gemmi/mmread_gz.hpp" // for read_structure_gz
#include "monlib_opt.h"

#define GEMMI_PROG crd
#include "options.h"

using namespace gemmi;

namespace {

enum OptionIndex {
  AutoCis=AfterMonLibOptions, AutoLink, AutoLigand,
  NoAliases, NoZeroOccRestr, NoHydrogens, KeepHydrogens
};

const option::Descriptor Usage[] = {
  { NoOp, 0, "", "", Arg::None,
    "Usage:"
    "\n " EXE_NAME " [options] INPUT_FILE OUTPUT_FILE"
    "\n\nPrepare intermediate Refmac files."
    "\nINPUT_FILE can be in PDB, mmCIF or mmJSON format."
    "\n\nOptions:" },
  CommonUsage[Help],
  CommonUsage[Version],
  CommonUsage[Verbose],
  MonLibUsage[0], // Monomers
  MonLibUsage[1], // Libin
  { AutoCis, 0, "", "auto-cis", Arg::YesNo,
    "  --auto-cis=Y|N  \tAssign cis/trans ignoring CISPEP record (default: Y)." },
  { AutoLink, 0, "", "auto-link", Arg::YesNo,
    "  --auto-link=Y|N  \tFind links not included in LINK/SSBOND (default: N)." },
  { AutoLigand, 0, "", "auto-ligand", Arg::YesNo,
    "  --auto-ligand=Y|N  \tIf ligand has no definition make ad-hoc restraints (N)." },
  { NoAliases, 0, "", "no-aliases", Arg::None,
    "  --no-aliases  \tIgnore _chem_comp_alias." },
  //{ NoZeroOccRestr, 0, "", "no-zero-occ", Arg::None,
  //  "  --no-zero-occ  \tNo restraints for zero-occupancy atoms." },
  { NoOp, 0, "", "", Arg::None,
    "\nHydrogen options (default: remove and add on riding positions):" },
  { NoHydrogens, 0, "H", "no-hydrogens", Arg::None,
    "  -H, --no-hydrogens  \tRemove (and do not add) hydrogens." },
  { KeepHydrogens, 0, "", "keep-hydrogens", Arg::None,
    "  --keep-hydrogens  \tPreserve hydrogens from the input file." },
  MonLibUsage[2], // details about Libin (--lib)
  { 0, 0, 0, 0, 0, 0 }
};

} // anonymous namespace

int GEMMI_MAIN(int argc, char **argv) {
  OptParser p(EXE_NAME);
  p.simple_parse(argc, argv, Usage);
  p.require_positional_args(2);
  p.check_exclusive_pair(KeepHydrogens, NoHydrogens);
  // check if monomer dir is specified (--monomers or $CLIBD_MON)
  MonArguments mon_args = get_monomer_args(p.options);
  std::string input = p.coordinate_input_file(0);
  std::string output = p.nonOption(1);
  int verbose = p.options[Verbose].count();

  try {
    if (verbose)
      fprintf(stderr, "Reading %s ...\n", input.c_str());
    cif::Document st_doc;
    Structure st = read_structure_gz(input, CoorFormat::Detect, &st_doc);
    setup_for_crd(st);

    if (st.models.empty()) {
      fprintf(stderr, "No models found in the input file.\n");
      return 1;
    }
    Model& model0 = st.models[0];
    std::vector<std::string> wanted = model0.get_all_residue_names();

    MonLib monlib;
    read_monomer_lib_and_user_files(monlib, wanted, mon_args, &st_doc);

    if (!wanted.empty()) {
      for (const std::string& name : wanted)
        fprintf(stderr, "WARNING: definition not found for %s.\n", name.c_str());
      if (!p.is_yes(AutoLigand, false))
        fail("Supply missing monomer definitions or use option --auto-ligand=Y");
      fprintf(stderr, "Note: Using ad-hoc restraints for missing monomers.\n"
                      "      Consider generating monomer CIFs with AceDRG or GRADE.\n");
    }

    if (p.options[NoAliases])
      for (auto& name_monomer : monlib.monomers)
        name_monomer.second.aliases.clear();

    bool use_cispeps = !p.is_yes(AutoCis, true);

    if (p.is_yes(AutoLink, false)) {
      size_t before = st.connections.size();
      add_automatic_links(model0, st, monlib);
      if (verbose)
        for (size_t i = before; i < st.connections.size(); ++i) {
          const Connection& conn = st.connections[i];
          fprintf(stderr, "Automatic link: %s - %s\n",
                  conn.partner1.str().c_str(), conn.partner2.str().c_str());
        }
    }

    if (verbose)
      fprintf(stderr, "Preparing topology, hydrogens, restraints...\n");
    bool reorder = true;
    bool ignore_unknown_links = false;
    HydrogenChange h_change;
    if (p.options[NoHydrogens])
      h_change = HydrogenChange::Remove;
    else if (p.options[KeepHydrogens])
      h_change = HydrogenChange::NoChange;
    else
      h_change = HydrogenChange::ReAddButWater;
    auto topo = prepare_topology(st, monlib, 0, h_change, reorder,
                                 &std::cerr, ignore_unknown_links, use_cispeps);
    if (!use_cispeps)
      topo->set_cispeps_in_structure(st);
    if (verbose)
      fprintf(stderr, "Preparing data for Refmac...\n");
    cif::Document crd = prepare_refmac_crd(st, *topo, monlib, h_change);
    // expand the starting comment
    cif::Item& first_item = crd.blocks.at(0).items.at(0);
    if (first_item.type == cif::ItemType::Comment) {
      std::string& comment = first_item.pair[1];
      comment += "\n# Command line: " EXE_NAME;
      for (int i = 1; i < argc; ++i)
        comment.append("  ").append(argv[i]);
    }
    if (verbose)
      fprintf(stderr, "Writing %s\n", output.c_str());
    Ofstream os(output, &std::cout);
    cif::WriteOptions cif_options;
    cif_options.compact = true;
    write_cif_to_stream(os.ref(), crd, cif_options);
  } catch (std::exception& e) {
    fprintf(stderr, "ERROR: %s\n", e.what());
    return 1;
  }
  return 0;
}