diff --git a/README.md b/README.md index 54e46d3b81..ff8eaaf0c6 100644 --- a/README.md +++ b/README.md @@ -196,7 +196,7 @@ Enhancements to entropy calculation in original Action\_Gist. Code for the CUDA version of the 'closest' Action. * Nathan D. Levinzon (University of Utah, UT) -Action\surftension. +Action\_surftension and TIMap. * Andrew Simmonett (National Institutes of Health) Code for the reciprocal part of the particle mesh Ewald calculation (electrostatic and Lennard-Jones). diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index dd657e6d8c..a57985356d 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -16345,7 +16345,7 @@ The following topology related commands are available: \begin_layout Standard \align center \begin_inset Tabular - + @@ -16801,6 +16801,26 @@ solvent Change which molecules are considered solvent. \end_layout +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +timap +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Map a residue onto a user template and write a sorted Amber OFF library; optional dual-topology TI pair. +\end_layout + \end_inset @@ -19337,6 +19357,828 @@ none is specified. \end_layout +\begin_layout Subsection +timap +\end_layout + +\begin_layout LyX-Code +timap +\begin_inset Index idx +range none +pageformat default +status collapsed + +\begin_layout Plain Layout +timap +\end_layout + +\end_inset + + [template ] [out ] [tiout ] +\begin_inset Newline newline +\end_inset + + [mapout ] [name ] +\end_layout + +\begin_layout LyX-Code + [maponly] [replace] [naorder] [aaorder] +\end_layout + +\begin_layout LyX-Code + [seed {auto|names|na|aa|none}] [anchor ] +\end_layout + +\begin_layout LyX-Code + [series ...] [outprefix ] [tioutprefix ] +\end_layout + +\begin_layout LyX-Code + [mapoutprefix ] [parentout ] +\end_layout + +\begin_deeper +\begin_layout Standard +Align atoms in topology +\series bold + +\series default + to a user-supplied template so shared atoms occupy the same indices. + By default this writes +\series bold + +\series default + as an Amber OFF +\series bold +.lib +\series default + in template atom order (no dummy atoms). + Dual-topology TI files are written only if +\series bold +tiout +\series default + is given. + This is intended for preparing residue libraries for thermodynamic integration (TI) of nucleic acids, amino acids, and small molecules, + where two (or more) end states must list common atoms in the same order. + For a series of analogs with no predetermined parent, +\series bold +series +\series default + auto-selects a minimum parent that shares the most with the rest, then maps every analog onto it. + Unlike +\series bold +atommap +\series default +, a one-to-one correspondence is not required: extra atoms on the target are insertions. + Partial maps are expected and useful. +\end_layout + +\begin_layout Standard +The template atom order +\series emph +is +\series default + the shared-atom order used for mapping. + Any topology can be a template (mol2, Amber OFF unit loaded with +\series bold +readdata +\series default +, and so on). + Example ModXNA fragments used by the tests live in +\series bold +test/Test_TIMap/ +\series default + (see PARENTS.txt). + The command does not load them. +\end_layout + +\begin_layout Description + Topology (or COORDS set) whose atoms will be reordered. +\end_layout + +\begin_layout Description +template +\begin_inset space ~ +\end_inset + + Parent topology whose +\series emph +current +\series default + atom order is the shared-atom order. + If omitted, +\series bold + +\series default + is mapped to itself; combine with +\series bold +naorder +\series default + to invent Amber's nucleic-acid walk with no parent file. +\end_layout + +\begin_layout Description +out +\begin_inset space ~ +\end_inset + + Write +\series bold + +\series default + as an Amber OFF library in template atom order. + Default (when +\series bold +maponly +\series default + is not given): +\series bold +.sorted.lib +\series default +. +\end_layout + +\begin_layout Description +mapout +\begin_inset space ~ +\end_inset + + Write the atom map (output index, target atom, template atom, insertion flag). +\end_layout + +\begin_layout Description +name +\begin_inset space ~ +\end_inset + + Store the remapped topology as a new topology set. +\end_layout + +\begin_layout Description +maponly Write the map (and +\series bold +tiout +\series default + files, if requested) but do not write the default aligned +\series bold +.lib +\series default + and do not store a remapped topology. +\end_layout + +\begin_layout Description +replace Overwrite +\series bold + +\series default + with the remapped atom order. +\end_layout + +\begin_layout Description +naorder Walk the template as P +\begin_inset Formula $\rightarrow$ +\end_inset + + O5' +\begin_inset Formula $\rightarrow$ +\end_inset + + C5' +\begin_inset Formula $\rightarrow$ +\end_inset + + C4' +\begin_inset Formula $\rightarrow$ +\end_inset + + O4' +\begin_inset Formula $\rightarrow$ +\end_inset + + C1' +\begin_inset Formula $\rightarrow$ +\end_inset + + base +\begin_inset Formula $\rightarrow$ +\end_inset + + C3' +\begin_inset Formula $\rightarrow$ +\end_inset + + C2' +\begin_inset Formula $\rightarrow$ +\end_inset + + O3' . + No parent is required: +\series bold + +\series default + is walked as itself. + If +\series bold +template +\series default + is also given, that parent is walked instead, then +\series bold + +\series default + is mapped onto it. +\end_layout + +\begin_layout Description +aaorder Walk +\series bold + +\series default + (or the template, if given) in ff19SB +\series bold +amino19.lib +\series default + order: N +\begin_inset Formula $\rightarrow$ +\end_inset + + H +\begin_inset Formula $\rightarrow$ +\end_inset + + CA +\begin_inset Formula $\rightarrow$ +\end_inset + + HA +\begin_inset Formula $\rightarrow$ +\end_inset + + side chain from CB (each heavy atom, then its hydrogens; smaller branch first) +\begin_inset Formula $\rightarrow$ +\end_inset + + C +\begin_inset Formula $\rightarrow$ +\end_inset + + O. + Proline is N +\begin_inset Formula $\rightarrow$ +\end_inset + + CD +\begin_inset Formula $\rightarrow$ +\end_inset + + ... +\begin_inset Formula $\rightarrow$ +\end_inset + + CB +\begin_inset Formula $\rightarrow$ +\end_inset + + CA +\begin_inset Formula $\rightarrow$ +\end_inset + + C +\begin_inset Formula $\rightarrow$ +\end_inset + + O. + No parent is required. + Use this to put a noncanonical amino-acid +\series bold +.lib +\series default + into the same shared-atom order as Amber ff19SB for TI. + The walk is built into the command; it does not read +\series bold +amino19.lib +\series default +. + A smaller amino-acid library already in the test tree is +\series bold +test/Test_Build/amino12.lib +\series default +. +\end_layout + +\begin_layout Description +seed How the correspondence is started. + +\series bold +auto +\series default + (default) uses nucleic-acid scaffold roles when they are detected, otherwise amino-acid N/CA/C/O/CB when a peptide motif is found, otherwise unique atom names. + +\series bold +names +\series default + uses unique names only. + +\series bold +na +\series default + forces nucleic-acid scaffold seeding. + +\series bold +aa +\series default + forces amino-acid N/CA/C/O/CB seeding. + +\series bold +none +\series default + performs no matching (not generally useful). + The map is then grown by neighbor element/degree signatures. + Partial maps are expected. +\end_layout + +\begin_layout Description +anchor +\begin_inset space ~ +\end_inset + + Name of the template atom before which completely unmatched leftovers are inserted (default O3'). + Atoms that are bonded to a mapped parent are inserted immediately after that parent instead. +\end_layout + +\begin_layout Description +tiout +\begin_inset space ~ +\end_inset + + Optional. + Write a dual-topology TI pair of equal size. + Lambda 0 is the template with charge-0 dummy atoms for target insertions; + lambda 1 is the target with charge-0 dummy atoms for unmatched template atoms. + Shared atoms occupy the same indices in both files. + Writes +\series bold +.0.mol2 +\series default + / +\series bold +.0.lib +\series default +, +\series bold +.1.mol2 +\series default + / +\series bold +.1.lib +\series default +, +\series bold +.scmask +\series default + (pmemd scmask1 / scmask2), and +\series bold +.atoms +\series default +. + Dummy atoms copy the partner's name and coordinates, have charge 0, mass 0, and Amber type DUM. + Does not suppress the default aligned +\series bold +.lib +\series default + unless +\series bold +maponly +\series default + is also given. + A remapped target topology is stored only if +\series bold +name +\series default + or +\series bold +replace +\series default + is also given. +\end_layout + +\begin_layout Description +series +\begin_inset space ~ +\end_inset + + ... List two or more analog topologies. + Chooses one parent that shares the most with the rest, then maps every listed analog onto that parent (star topology for TI). + Do not combine with +\series bold +template +\series default +. + Parent selection: maximize the sum of mapped atoms when every other member is matched onto the candidate; on ties prefer the fewest atoms, then a leadopt-style similarity sum +\begin_inset Formula $\sum_{j}\exp(-0.1\cdot(\mathrm{insertions}_{j}+\mathrm{unmapped}_{j}))$ +\end_inset + +. + That score is adapted from +\series bold +util/leadopt/similarity.py +\series default + ( +\series emph +exp_delta +\series default + / +\series emph +by_heavy_atom_count +\series default +; Summa Lab / Jonathan Redmann & Christopher Summa, University of New Orleans). + The common substructure is approximated by +\series bold +TIMatch +\series default +'s partial map rather than an external MCS engine ( +\series bold +util/leadopt/mcs.py +\series default +, +\series bold +rule.py +\series default +, +\series bold +graph.py +\series default +). + The star-parent use case is inspired by +\series bold +util/leadopt/GraphGenerator4.py +\series default + (Redmann & Summa) without importing that TI-path planner. + The +\series bold +leadopt/ +\series default + tree is reference-only and is not linked into the cpptraj build. +\end_layout + +\begin_layout Description +outprefix +\begin_inset space ~ +\end_inset + + With +\series bold +series +\series default +, write each aligned OFF library as +\series bold +.sorted.lib +\series default + (default: +\series bold +.sorted.lib +\series default +). +\end_layout + +\begin_layout Description +tioutprefix +\begin_inset space ~ +\end_inset + + With +\series bold +series +\series default +, write dual-topology TI files for every non-parent member as +\series bold +.{0,1}.{mol2,lib} +\series default +, +\series bold +.scmask +\series default +, and +\series bold +.atoms +\series default +. +\end_layout + +\begin_layout Description +mapoutprefix +\begin_inset space ~ +\end_inset + + With +\series bold +series +\series default +, write each atom map as +\series bold +.map +\series default +. +\end_layout + +\begin_layout Description +parentout +\begin_inset space ~ +\end_inset + + With +\series bold +series +\series default +, write the selected parent name and per-candidate scores. +\end_layout + +\begin_layout Standard +\series bold +Generating your own templates +\series default +. + If +\series bold +template +\series default + is given, that parent's current atom order is used. + If it is omitted, +\series bold +naorder +\series default + or +\series bold +aaorder +\series default + invents Amber's canonical walk on +\series bold + +\series default + itself. +\end_layout + +\begin_layout LyX-Code +# 1. Parent given: analog is sorted to that residue's atom order. +\end_layout + +\begin_layout LyX-Code +# 2. No parent: naorder / aaorder invents Amber's canonical walk on the analog. +\end_layout + +\begin_layout LyX-Code +# 3. Optional tiout writes dual-topology TI files. +\end_layout + +\begin_layout Standard +Canonical NA order, no parent (sort this residue the way Amber walks a nucleotide). +\end_layout + +\begin_layout LyX-Code +> readdata FLE.lib name FLE +\end_layout + +\begin_layout LyX-Code +> timap FLE[FLE] naorder out FLE.lib +\end_layout + +\begin_layout LyX-Code +> go +\end_layout + +\begin_layout Standard +Nucleic acids onto a parent (example: FLE = unmodified rA onto ERN = unmodified dA). + Shared atoms keep the dA indices; the 2'-OH is an insertion next to C2'; dA H2'' is unmatched on the template. + Default: write FLE in ERN atom order as an OFF library. +\end_layout + +\begin_layout LyX-Code +> readdata ERN.lib name ERN +\end_layout + +\begin_layout LyX-Code +> readdata FLE.lib name FLE +\end_layout + +\begin_layout LyX-Code +> timap FLE[FLE] template ERN[ERN] out FLE.lib +\end_layout + +\begin_layout LyX-Code +> go +\end_layout + +\begin_layout Standard +Add +\series bold +tiout +\series default + to also write dual-topology libraries (34 atoms each) and pmemd scmask lines. +\end_layout + +\begin_layout LyX-Code +> timap FLE[FLE] template ERN[ERN] tiout fle_ern +\end_layout + +\begin_layout LyX-Code +> go +\end_layout + +\begin_layout Standard +Small molecules (example: phenol onto benzene). + The ring maps by unique names; OH inserts after C1; benzene H1 is unmatched. +\end_layout + +\begin_layout LyX-Code +> parm benzene.mol2 name benzene +\end_layout + +\begin_layout LyX-Code +> parm phenol.mol2 name phenol +\end_layout + +\begin_layout LyX-Code +> timap phenol template benzene out phenol.lib +\end_layout + +\begin_layout LyX-Code +> go +\end_layout + +\begin_layout Standard +Amino acids (example: selenocysteine onto cysteine). + Use Amber OFF so selenium is +\series bold +elmnt 34 +\series default +. + Mol2 atom name SE is read as sulfur. + Backbone names map; SE is an insertion at CB; cysteine SG is unmatched. +\end_layout + +\begin_layout LyX-Code +> readdata cys.lib name CYS +\end_layout + +\begin_layout LyX-Code +> readdata sec.lib name SEC +\end_layout + +\begin_layout LyX-Code +> timap SEC[SEC] template CYS[CYS] out SEC.lib +\end_layout + +\begin_layout LyX-Code +> go +\end_layout + +\begin_layout Standard +Canonical amino-acid order, no parent (sort a noncanonical +\series bold +.lib +\series default + the way ff19SB walks Cys/Ala/...). +\end_layout + +\begin_layout LyX-Code +> readdata ncaa.lib name ncaa +\end_layout + +\begin_layout LyX-Code +> timap ncaa[ncaa] aaorder out ncaa.lib +\end_layout + +\begin_layout LyX-Code +> go +\end_layout + +\begin_layout Standard +Series of analogs (auto-select parent, then map each onto it for TI). + Example: phenol and benzene; benzene is chosen (same mapped total, fewer atoms). +\end_layout + +\begin_layout LyX-Code +> parm benzene.mol2 name benzene +\end_layout + +\begin_layout LyX-Code +> parm phenol.mol2 name phenol +\end_layout + +\begin_layout LyX-Code +> timap series phenol benzene outprefix series_ tioutprefix series_ti_ parentout series_parent.dat +\end_layout + +\begin_layout LyX-Code +> go +\end_layout + +\begin_layout Subsubsection* +How the nucleic-acid walk was chosen +\end_layout + +\begin_layout Standard +ModXNA ( +https://github.com/modna-ff +) is the baseline chemistry for +\series bold +naorder +\series default +. + It parameterizes nucleic acids as fragments (nucleobase, sugar, backbone) rather than as whole residues. + Combining those fragments into one +\series bold +.lib +\series default + does not put shared atoms in a stable order. + That breaks thermodynamic integration (common atoms must occupy the same indices) and makes residues harder for +\series bold +nastruct +\series default + to treat as nucleotides, because +\series bold +nastruct +\series default + identifies sugars and bases from connectivity and Amber-like names (C1', N9/N1, and so on), not from file order. +\end_layout + +\begin_layout Standard +The walk is therefore not copied from any one ModXNA file. + The fragment library was grouped by chemical class (backbone with or without phosphate, purine vs pyrimidine base, 2' pattern of the sugar). + One consensus parent was kept for each class; those 11 example fragments are in +\series bold +test/Test_TIMap/ +\series default + (see PARENTS.txt). + They are test data, not installed with cpptraj. + Each parent's atom list is the recommended shared-atom order for mapping other members of that class. + Combined nucleotides are a separate case: they are mapped to a user residue (for example unmodified rA onto dA), not to those fragments. +\end_layout + +\begin_layout Standard +The +\series bold +naorder +\series default + walk is a 5' +\begin_inset Formula $\rightarrow$ +\end_inset + + 3' sugar-phosphate traversal chosen so that (1) the shared scaffold is a fixed prefix independent of how the file was written, (2) chemical insertions sit next to their attachment atom (nucleobase analogs after the glycosidic nitrogen +\begin_inset Formula $\chi$ +\end_inset + +, 2'-OH / OMe / F after C2'), and (3) the 3' linker O3' remains last among scaffold atoms so unmatched leftovers can be placed just before it. + Hydrogens follow their heavy atom, sorted by name, matching Amber convention. +\end_layout + +\begin_layout LyX-Code +P -> OP (name-sorted) -> O5' -> C5' -> C4' -> O4' -> C1' -> +\end_layout + +\begin_layout LyX-Code + base (BFS from chi, never recrossing C1') -> C3' -> C2' + 2' substituents -> O3' +\end_layout + +\begin_layout Standard +That is the order +\series bold +naorder +\series default + applies when no +\series bold +template +\series default + is given (the residue is walked as itself). + When a ModXNA parent +\series emph +is +\series default + given, +\series bold +timap +\series default + uses that parent's current atom order instead; those files were already written in the class order above. + The amino-acid walk ( +\series bold +aaorder +\series default +) plays the same role using ff19SB +\series bold +amino19.lib +\series default + as the baseline. +\end_layout + +\begin_layout Standard +Alias: +\series bold +timatch +\series default +. +\end_layout + +\end_deeper + \begin_layout Subsection updateparameters \end_layout @@ -22708,7 +23550,7 @@ Mod \begin_layout Standard \align center \begin_inset Tabular - + diff --git a/src/Command.cpp b/src/Command.cpp index e15a1edee1..bddaf600ab 100644 --- a/src/Command.cpp +++ b/src/Command.cpp @@ -73,6 +73,7 @@ #include "Exec_ParmStrip.h" #include "Exec_ParmWrite.h" #include "Exec_ScaleDihedralK.h" +#include "Exec_TIMap.h" #include "Exec_Top.h" #include "Exec_UpdateParameters.h" #include "Exec_Desc.h" @@ -328,6 +329,7 @@ void Command::Init() { Command::AddCmd( new Exec_ParmWrite(), Cmd::EXE, 1, "parmwrite" ); Command::AddCmd( new Exec_ResInfo(), Cmd::EXE, 1, "resinfo" ); Command::AddCmd( new Exec_ScaleDihedralK(),Cmd::EXE, 1, "scaledihedralk" ); + Command::AddCmd( new Exec_TIMap(), Cmd::EXE, 2, "timap", "timatch" ); Command::AddCmd( new Exec_UBInfo(), Cmd::EXE, 2, "ubinfo", "printub" ); Command::AddCmd( new Exec_UpdateParameters(), Cmd::EXE, 1, "updateparameters"); // ACTION diff --git a/src/Exec_TIMap.cpp b/src/Exec_TIMap.cpp new file mode 100644 index 0000000000..0e3def7160 --- /dev/null +++ b/src/Exec_TIMap.cpp @@ -0,0 +1,1314 @@ +#include "Exec_TIMap.h" +#include "TIMatch.h" +#include "CpptrajStdio.h" +#include "CpptrajFile.h" +#include "CpptrajState.h" +#include "DataSet_Topology.h" +#include "DataSet_Coords.h" +#include "DataFile.h" +#include "ParmFile.h" +#include "FileName.h" +#include "StringRoutines.h" +#include "ArgList.h" +#include "Frame.h" +#include "CoordinateInfo.h" +#include "Trajout_Single.h" +#include "Trajin_Single.h" +#include "TrajectoryFile.h" +#include "Residue.h" +#include "Atom.h" +#include "Topology.h" +#include +#include +#include +#include + +/* Source encoding: UTF-8. Unicode (O5′, 2′, λ) is comments-only. + * + * Workflow: + * 1. Match target onto template (partial map). + * 2. Write the remapped target (out / fmt: Amber OFF .lib, mol2, or pdb). + * 3. If tiout is given, also build a dual-topology pair of the + * same size: unique atoms are real in one end state and charge-0 / + * mass-0 / type-DUM dummies in the other. + * 4. series mode: choose one parent among many analogs, then map every + * other analog onto that parent. + * + * Series scoring uses exp(-0.1 * (insertions + unmatched template atoms)), + * the similarity form from leadopt (Redmann & Summa, University of New + * Orleans). The common substructure is TIMatch's partial map, not an + * external MCS engine. + * \author Nathan D. Levinzon + */ + +// Exec_TIMap::Help() +void Exec_TIMap::Help() const { + mprintf("\t [template ] [naorder|aaorder]\n" + "\t[out ] [fmt {lib|mol2|pdb}] [tiout ] [mapout ] [maponly]\n" + "\t[name ] [replace] [seed {auto|names|na|aa|none}] [anchor ]\n" + "\t[series ... [outprefix

] [tioutprefix

] [mapoutprefix

]\n" + "\t [parentout ]]\n" + " Reorder so atoms shared with a template have the same indices (TI).\n" + " Partial maps are expected (unlike atommap). Alias: timatch.\n" + " Type 'help timap extended' for keywords, output files, and examples.\n"); +} + +// Exec_TIMap::Help() +void Exec_TIMap::Help(ArgList& argIn) const { + if (!argIn.hasKey("extended")) { + Help(); + return; + } + mprintf("\t{ [template ] [naorder|aaorder] |\n" + "\t series ... [outprefix ] [tioutprefix ]\n" + "\t [mapoutprefix ] [parentout ] }\n" + "\t[out ] [fmt {lib|mol2|pdb}] [tiout ] [mapout ] [maponly]\n" + "\t[name ] [replace]\n" + "\t[seed {auto|names|na|aa|none}] [anchor ]\n" + " Reorder atoms so shared sites share indices for TI (nucleotides,\n" + " amino acids, small molecules). Partial maps are expected (unlike atommap).\n" + "\n" + " Atom order (pick one; combine template with naorder/aaorder to walk the\n" + " parent first when its file order is messy):\n" + " template parent's current order (Amber .lib, mol2, ...).\n" + " naorder Amber NA walk of (or of template if given):\n" + " P -> OP -> O5' -> C5' -> C4' -> O4' -> C1' -> base ->\n" + " C3' -> C2' -> O3'.\n" + " aaorder ff19SB amino19.lib walk of (or of template):\n" + " N -> H -> CA -> HA -> side chain from CB -> C -> O\n" + " (proline: N -> CD -> ... -> CB -> CA -> C -> O).\n" + " Does not read amino19.lib; the walk is built in.\n" + " series ... auto-pick one parent among the listed topologies\n" + " (max mapped-atom sum; ties: fewest atoms, then\n" + " exp(-0.1*(insertions+unmatched)) sum),\n" + " then map every member onto it. Not with template.\n" + "\n" + " Output:\n" + " Default: remapped structure (out ; else .sorted.lib).\n" + " fmt {lib|mol2|pdb} choose writer (default lib). If out has extension\n" + " .mol2 / .pdb / .lib, that format is used instead.\n" + " maponly skip the aligned structure / in-memory parm; still write\n" + " mapout/tiout.\n" + " mapout atom correspondence file.\n" + " name / replace store remapped topology, or overwrite .\n" + " tiout dual-topology TI pair (same NATOM; DUM q=0 m=0):\n" + " .0.mol2/.lib lambda=0 template + dummy insertions\n" + " .1.mol2/.lib lambda=1 target + dummy unmatched atoms\n" + " .scmask pmemd scmask1 / scmask2\n" + " .atoms per-slot kind, names, charges\n" + " series also: .sorted.{lib|mol2|pdb},\n" + " .*, .map,\n" + " parentout score table.\n" + "\n" + " Matching: seed auto|names|na|aa|none (default auto). Insertions bonded\n" + " to a mapped parent follow that parent; other leftovers go before\n" + " anchor (default O3').\n" + "\n" + " Inputs: Amber OFF .lib (readdata -> Name[Unit]), or mol2 / pdb / Amber\n" + " topology via parm — or pass a file path and timap loads it. Formats may\n" + " be mixed; everything is matched as Topology. Supply your own template.\n" + " Immediate command; add go if the input has trajin. Alias: timatch.\n" + "\n" + " Examples:\n" + " > readdata parent.lib name parent\n" + " > readdata analog.lib name analog\n" + " > timap analog[analog] template parent[parent] out analog.lib\n" + " > timap analog[analog] naorder out analog.lib\n" + " > timap ncaa[ncaa] aaorder out ncaa.lib\n" + " > timap analog[analog] template parent[parent] tiout analog_ti\n" + " > parm a.mol2 name A\n" + " > parm b.pdb name B\n" + " > timap series A B outprefix s_ tioutprefix s_ti_ parentout parent.dat\n" + " > timap phenol.mol2 template benzene.pdb out phenol.lib\n" + " > timap phenol.mol2 template cys.lib[CYS] mapout mix.map maponly\n"); +} + +/** Look up a topology by the name the user typed (already-loaded sets only). + * + * Order of attempts: + * 1. A TOPOLOGY dataset (parm / parmwrite sets: mol2, pdb, Amber, ...). + * 2. A COORDS dataset — Amber OFF units loaded with readdata appear as + * LibName[UnitName], e.g. FLE[FLE] or CYS[CYS]. + * 3. A purely numeric string, treated as a parm index. + * + * \param dsOut If non-null, receives the dataset pointer (TOPOLOGY or COORDS). + */ +Topology* Exec_TIMap::FindNamedTop(DataSetList& dsl, std::string const& name, DataSet** dsOut) const +{ + if (dsOut != 0) *dsOut = 0; + if (name.empty()) return 0; + DataSet* ds = dsl.FindSetOfType(name, DataSet::TOPOLOGY); + if (ds != 0) { + if (dsOut != 0) *dsOut = ds; + return ((DataSet_Topology*)ds)->TopPtr(); + } + ds = dsl.FindSetOfGroup(name, DataSet::COORDINATES); + if (ds != 0) { + if (dsOut != 0) *dsOut = ds; + return ((DataSet_Coords*)ds)->TopPtr(); + } + bool isInt = !name.empty(); + for (size_t i = 0; i < name.size(); i++) { + if (name[i] < '0' || name[i] > '9') { isInt = false; break; } + } + if (!isInt) return 0; + ArgList tmp(name); + return dsl.GetTopByIndex(tmp); +} + +/** Split optional Amber-OFF unit bracket: path[Unit] → path, Unit. */ +void Exec_TIMap::ParseTopSpec(std::string const& spec, std::string& primary, std::string& unit) const +{ + primary = spec; + unit.clear(); + if (spec.empty()) return; + size_t rb = spec.rfind(']'); + size_t lb = spec.rfind('['); + if (lb != std::string::npos && rb != std::string::npos && + rb == spec.size() - 1 && rb > lb) + { + primary = spec.substr(0, lb); + unit = spec.substr(lb + 1, rb - lb - 1); + } +} + +/** Basename without extension (phenol.mol2 → phenol). */ +std::string Exec_TIMap::FileStem(FileName const& fn) const +{ + std::string b = fn.Base(); + std::string e = fn.Ext(); + if (!e.empty() && b.size() > e.size() && + b.compare(b.size() - e.size(), e.size(), e) == 0) + return b.substr(0, b.size() - e.size()); + return b; +} + +/** + * Resolve a user topology specifier to Topology. + * Accepts an already-loaded set name, or a mol2/pdb/Amber/.lib file path + * (auto-loaded via ParmFile or DataFile Amber OFF). All formats become the + * same internal Topology for TIMatch; mix-and-match is allowed. + */ +Topology* Exec_TIMap::ResolveTop(CpptrajState& State, std::string const& spec, DataSet** dsOut) const +{ + if (dsOut != 0) *dsOut = 0; + if (spec.empty()) return 0; + + Topology* top = FindNamedTop(State.DSL(), spec, dsOut); + if (top != 0) return top; + + std::string primary, unit; + ParseTopSpec(spec, primary, unit); + if (primary != spec) { + top = FindNamedTop(State.DSL(), primary, dsOut); + if (top != 0) return top; + if (!unit.empty()) { + std::string withUnit = primary + "[" + unit + "]"; + top = FindNamedTop(State.DSL(), withUnit, dsOut); + if (top != 0) return top; + } + } + + bool maybeFile = File::Exists(primary) || + primary.find('.') != std::string::npos || + primary.find('/') != std::string::npos || + primary.find('\\') != std::string::npos; + if (!maybeFile || !File::Exists(primary)) + return 0; + + FileName fn(primary); + std::string stem = FileStem(fn); + if (stem.empty()) stem = fn.Base(); + + // Topology formats: mol2, pdb, Amber parm, PSF, ... + if (ParmFile::DetectFormat(fn) != ParmFile::UNKNOWN_PARM) { + ArgList nameArg(std::string("name ") + stem); + if (State.AddTopology(primary, nameArg)) { + mprinterr("Error: timap: could not load topology file '%s'\n", primary.c_str()); + return 0; + } + mprintf("\tLoaded topology file '%s' as '%s'\n", primary.c_str(), stem.c_str()); + return FindNamedTop(State.DSL(), stem, dsOut); + } + + // Amber OFF .lib (and other DataFile formats that yield COORDS with Topology) + DataFile dfile; + ArgList nameArg(std::string("name ") + stem); + if (dfile.ReadDataIn(fn, nameArg, State.DSL())) { + mprinterr("Error: timap: could not load data/topology file '%s'\n", primary.c_str()); + return 0; + } + std::string lookup = unit.empty() ? stem : (stem + "[" + unit + "]"); + top = FindNamedTop(State.DSL(), lookup, dsOut); + if (top == 0 && !unit.empty()) + top = FindNamedTop(State.DSL(), stem, dsOut); + if (top == 0) { + mprinterr("Error: timap: loaded '%s' but found no topology set named '%s'\n", + primary.c_str(), lookup.c_str()); + return 0; + } + mprintf("\tLoaded file '%s' as '%s'\n", primary.c_str(), lookup.c_str()); + return top; +} + +/** Fill frm with coordinates for top. Prefer a COORDS frame; else the original file. */ +bool Exec_TIMap::LoadCoords(Topology const& top, DataSet* ds, Frame& frm) const +{ + frm.SetupFrame(top.Natom()); + for (int i = 0; i < frm.Natom() * 3; i++) + frm[i] = 0.0; + if (ds != 0 && ds->Group() == DataSet::COORDINATES && ds->Size() > 0) { + DataSet_Coords* crd = (DataSet_Coords*)ds; + Frame tmp = crd->AllocateFrame(); + crd->GetFrame(0, tmp); + int ncopy = std::min(tmp.Natom(), frm.Natom()); + for (int i = 0; i < ncopy * 3; i++) + frm[i] = tmp[i]; + return ncopy == frm.Natom(); + } + if (top.OriginalFilename().empty()) + return frm.Natom() == 0; + Trajin_Single traj; + ArgList empty; + Topology* tptr = const_cast(&top); + if (traj.SetupTrajRead(top.OriginalFilename(), empty, tptr)) + return false; + if (traj.BeginTraj()) + return false; + int got = traj.GetNextFrame(frm); + traj.EndTraj(); + return got == 1; +} + +void Exec_TIMap::CopyXyz(Frame const& src, int oldat, Frame& dst, int newat) const +{ + if (oldat < 0 || newat < 0) return; + if (oldat >= src.Natom() || newat >= dst.Natom()) return; + dst[newat * 3 ] = src[oldat * 3 ]; + dst[newat * 3 + 1] = src[oldat * 3 + 1]; + dst[newat * 3 + 2] = src[oldat * 3 + 2]; +} + +/** Permute src coordinates into dst using Map[new] = old. */ +void Exec_TIMap::ApplyOrderToFrame(Frame const& src, std::vector const& order, Frame& dst) const +{ + dst.SetupFrame((int)order.size()); + for (int i = 0; i < (int)order.size() * 3; i++) + dst[i] = 0.0; + for (int i = 0; i < (int)order.size(); i++) + CopyXyz(src, order[i], dst, i); +} + +/** Amber OFF unit names: letters, digits, underscore only. */ +std::string Exec_TIMap::OffUnitName(std::string in, char fallback) const +{ + std::string out; + for (size_t i = 0; i < in.size(); i++) { + char c = in[i]; + if (std::isalnum((unsigned char)c) || c == '_') + out += c; + } + if (out.empty()) { + out = "L"; + out += fallback; + } + if (out.size() > 20) out.resize(20); + return out; +} + +NameType Exec_TIMap::ResNameOf(Topology const& top, const char* fallback) const +{ + if (top.Nres() > 0) + return top.Res(0).Name(); + return NameType(fallback); +} + +/** Build one TI end-state topology+frame from the dual-slot list. + * lambda0: template atoms real, target insertions dummy (charge 0, mass 0, type DUM). + * lambda1: target atoms real, unmatched template atoms dummy. + */ +int Exec_TIMap::BuildTiUnit(bool lambda0, + Topology const& tgt, Topology const& tpl, + Frame const& tgtX, Frame const& tplX, + std::vector const& dual, + Topology& outTop, Frame& outFrm, + NameType const& resname) const +{ + typedef TIMatch::Result::DualSlot Slot; + outTop = Topology(); + Residue res(resname, 1, ' ', ""); + int nd = (int)dual.size(); + outFrm.SetupFrame(nd); + for (int i = 0; i < nd * 3; i++) + outFrm[i] = 0.0; + std::vector tgtDual(tgt.Natom(), -1); + std::vector tplDual(tpl.Natom(), -1); + for (int i = 0; i < nd; i++) { + Slot const& s = dual[i]; + if (s.tgt_ >= 0) tgtDual[s.tgt_] = i; + if (s.tpl_ >= 0) tplDual[s.tpl_] = i; + Atom src; + bool dummy = false; + int xyzFrom = -1; + Frame const* xyzSrc = 0; + if (lambda0) { + if (s.tpl_ >= 0) { + src = tpl[s.tpl_]; + xyzSrc = &tplX; + xyzFrom = s.tpl_; + } else { + src = tgt[s.tgt_]; + xyzSrc = &tgtX; + xyzFrom = s.tgt_; + dummy = true; + } + } else { + if (s.tgt_ >= 0) { + src = tgt[s.tgt_]; + xyzSrc = &tgtX; + xyzFrom = s.tgt_; + } else { + src = tpl[s.tpl_]; + xyzSrc = &tplX; + xyzFrom = s.tpl_; + dummy = true; + } + } + src.ClearBonds(); + if (dummy) { + src.SetCharge(0.0); + src.SetMass(0.0); + src.SetTypeName("DUM"); + } + outTop.AddTopAtom(src, res); + if (xyzSrc != 0) + CopyXyz(*xyzSrc, xyzFrom, outFrm, i); + } + + if (lambda0) { + for (int r = 0; r < tpl.Natom(); r++) { + if (tplDual[r] < 0) continue; + Atom const& at = tpl[r]; + for (Atom::bond_iterator b = at.bondbegin(); b != at.bondend(); ++b) { + if (*b > r && tplDual[*b] >= 0) + outTop.AddBond(tplDual[r], tplDual[*b], -1); + } + } + for (int i = 0; i < nd; i++) { + if (!dual[i].IsTgtOnly()) continue; + int n = dual[i].tgt_; + int attach = -1; + Atom const& at = tgt[n]; + for (Atom::bond_iterator b = at.bondbegin(); b != at.bondend(); ++b) { + int d = tgtDual[*b]; + if (d < 0) continue; + if (attach < 0 || dual[d].IsShared()) + attach = d; + } + if (attach >= 0) + outTop.AddBond(i, attach, -1); + else + mprintf("Warning: timap: dummy %s has no bonded neighbor in the dual layout.\n", + outTop[i].c_str()); + } + } else { + for (int n = 0; n < tgt.Natom(); n++) { + if (tgtDual[n] < 0) continue; + Atom const& at = tgt[n]; + for (Atom::bond_iterator b = at.bondbegin(); b != at.bondend(); ++b) { + if (*b > n && tgtDual[*b] >= 0) + outTop.AddBond(tgtDual[n], tgtDual[*b], -1); + } + } + for (int i = 0; i < nd; i++) { + if (!dual[i].IsTplOnly()) continue; + int r = dual[i].tpl_; + int attach = -1; + Atom const& at = tpl[r]; + for (Atom::bond_iterator b = at.bondbegin(); b != at.bondend(); ++b) { + int d = tplDual[*b]; + if (d < 0) continue; + if (attach < 0 || dual[d].IsShared()) + attach = d; + } + if (attach >= 0) + outTop.AddBond(i, attach, -1); + else + mprintf("Warning: timap: dummy %s has no bonded neighbor in the dual layout.\n", + outTop[i].c_str()); + } + } + outTop.CommonSetup(); + return 0; +} + +int Exec_TIMap::WriteMol2File(std::string const& fname, Topology& top, Frame const& frm, + DataSetList const& dsl) const +{ + Trajout_Single out; + ArgList empty; + if (out.PrepareTrajWrite(fname, empty, dsl, &top, CoordinateInfo(), 1, + TrajectoryFile::MOL2FILE)) + { + mprinterr("Error: timap: could not set up mol2 '%s'\n", fname.c_str()); + return 1; + } + if (out.WriteSingle(0, frm)) { + mprinterr("Error: timap: writing mol2 '%s'\n", fname.c_str()); + return 1; + } + out.EndTraj(); + mprintf("\tWrote '%s' (%i atoms).\n", fname.c_str(), top.Natom()); + return 0; +} + +/** Write PDB via Traj_PDBfile; include CONECT so bonds round-trip for rematching. */ +int Exec_TIMap::WritePdbFile(std::string const& fname, Topology& top, Frame const& frm, + DataSetList const& dsl) const +{ + Trajout_Single out; + ArgList pdbArgs("conect"); + if (out.PrepareTrajWrite(fname, pdbArgs, dsl, &top, CoordinateInfo(), 1, + TrajectoryFile::PDBFILE)) + { + mprinterr("Error: timap: could not set up pdb '%s'\n", fname.c_str()); + return 1; + } + if (out.WriteSingle(0, frm)) { + mprinterr("Error: timap: writing pdb '%s'\n", fname.c_str()); + return 1; + } + out.EndTraj(); + mprintf("\tWrote '%s' (%i atoms).\n", fname.c_str(), top.Natom()); + return 0; +} + +const char* Exec_TIMap::AlignedOutExt(AlignedOutFmt f) const +{ + switch (f) { + case AOUT_MOL2: return ".sorted.mol2"; + case AOUT_PDB: return ".sorted.pdb"; + case AOUT_LIB: + default: return ".sorted.lib"; + } +} + +const char* Exec_TIMap::AlignedOutLabel(AlignedOutFmt f) const +{ + switch (f) { + case AOUT_MOL2: return "mol2"; + case AOUT_PDB: return "pdb"; + case AOUT_LIB: + default: return "Amber OFF .lib"; + } +} + +int Exec_TIMap::ParseAlignedOutFmt(std::string const& s, AlignedOutFmt& out) const +{ + std::string l = ToLower(s); + if (l == "lib" || l == "off" || l == "amberlib") { out = AOUT_LIB; return 0; } + if (l == "mol2") { out = AOUT_MOL2; return 0; } + if (l == "pdb") { out = AOUT_PDB; return 0; } + return 1; +} + +/** Prefer extension of out ; otherwise keep def (from fmt keyword). */ +Exec_TIMap::AlignedOutFmt Exec_TIMap::FmtFromFilename(std::string const& fname, AlignedOutFmt def) const +{ + FileName fn(fname); + std::string e = ToLower(fn.Ext()); + if (e == ".mol2") return AOUT_MOL2; + if (e == ".pdb") return AOUT_PDB; + if (e == ".lib" || e == ".off") return AOUT_LIB; + return def; +} + +std::string Exec_TIMap::Q(std::string const& s) const +{ + return std::string("\"") + s + "\""; +} + +std::string Exec_TIMap::AtomTypeStr(Atom const& a) const +{ + std::string t = a.Type().Truncated(); + if (t.empty() && a.ElementName() != 0) + t = a.ElementName(); + return t; +} + +/** 1-based index of the first atom with this name, or 0 if missing. */ +int Exec_TIMap::AtomNum1(Topology const& top, const char* n) const +{ + for (int i = 0; i < top.Natom(); i++) { + if (top[i].Name() == n) + return i + 1; + } + return 0; +} + +/** Head/tail connect atoms for LEaP (P/O3' or N/C). */ +void Exec_TIMap::LeapConnect(Topology const& top, int& c1, int& c2, const char*& restype) const +{ + c1 = AtomNum1(top, "P"); + c2 = AtomNum1(top, "O3'"); + if (c1 > 0 && c2 > 0) { + restype = "n"; + return; + } + c1 = AtomNum1(top, "N"); + c2 = AtomNum1(top, "C"); + if (c1 > 0 && c2 > 0) { + restype = "p"; + return; + } + c1 = 0; + c2 = 0; + restype = "?"; +} + +/** Minimal Amber OFF writer for a single-residue TI unit. */ +int Exec_TIMap::WriteAmberLib(std::string const& fname, std::string const& unit, + Topology const& top, Frame const& frm) const +{ + CpptrajFile out; + if (out.OpenWrite(fname)) { + mprinterr("Error: Could not open '%s' for Amber OFF write.\n", fname.c_str()); + return 1; + } + int n = top.Natom(); + std::string res = (top.Nres() > 0) ? top.Res(0).Name().Truncated() : unit; + int c1, c2; + const char* restype; + LeapConnect(top, c1, c2, restype); + out.Printf("!!index array str\n \"%s\"\n", unit.c_str()); + out.Printf("!entry.%s.unit.atoms table str name str type int typex " + "int resx int flags int seq int elmnt dbl chg\n", unit.c_str()); + for (int i = 0; i < n; i++) { + Atom const& a = top[i]; + out.Printf(" %s %s %5d %3d %7d %3d %3d %10.6f\n", + Q(a.Name().Truncated()).c_str(), + Q(AtomTypeStr(a)).c_str(), + 0, 1, 131072, i + 1, a.AtomicNumber(), a.Charge()); + } + out.Printf("!entry.%s.unit.atomspertinfo table str pname str ptype " + "int ptypex int pelmnt dbl pchg\n", unit.c_str()); + for (int i = 0; i < n; i++) + out.Printf(" %s %s 0 -1 0.0\n", + Q(top[i].Name().Truncated()).c_str(), + Q(AtomTypeStr(top[i])).c_str()); + out.Printf("!entry.%s.unit.boundbox array dbl\n -1.000000\n 0.0\n 0.0\n 0.0\n 0.0\n", + unit.c_str()); + out.Printf("!entry.%s.unit.childsequence single int\n 2\n", unit.c_str()); + out.Printf("!entry.%s.unit.connect array int\n %d\n %d\n", unit.c_str(), c1, c2); + out.Printf("!entry.%s.unit.connectivity table int atom1x int atom2x int flags\n", + unit.c_str()); + for (int i = 0; i < n; i++) { + Atom const& a = top[i]; + for (Atom::bond_iterator b = a.bondbegin(); b != a.bondend(); ++b) { + if (*b > i) + out.Printf(" %d %d 1\n", i + 1, (*b) + 1); + } + } + out.Printf("!entry.%s.unit.hierarchy table str abovetype int abovex " + "str belowtype int belowx\n", unit.c_str()); + out.Printf(" \"U\" 0 \"R\" 1\n"); + for (int i = 0; i < n; i++) + out.Printf(" \"R\" 1 \"A\" %d\n", i + 1); + out.Printf("!entry.%s.unit.name single str\n %s\n", unit.c_str(), Q(unit).c_str()); + out.Printf("!entry.%s.unit.positions table dbl x dbl y dbl z\n", unit.c_str()); + for (int i = 0; i < n; i++) + out.Printf(" %.6f %.6f %.6f\n", frm[i * 3], frm[i * 3 + 1], frm[i * 3 + 2]); + out.Printf("!entry.%s.unit.residueconnect table int c1x int c2x int c3x " + "int c4x int c5x int c6x\n %d %d 0 0 0 0\n", unit.c_str(), c1, c2); + out.Printf("!entry.%s.unit.residues table str name int seq int childseq " + "int startatomx str restype int imagingx\n", unit.c_str()); + out.Printf(" %s 1 %d 1 \"%s\" 0\n", Q(res).c_str(), n + 1, restype); + out.Printf("!entry.%s.unit.residuesPdbSequenceNumber array int\n 0\n", unit.c_str()); + out.Printf("!entry.%s.unit.solventcap array dbl\n -1.000000\n 0.0\n 0.0\n 0.0\n 0.0\n", + unit.c_str()); + out.Printf("!entry.%s.unit.velocities table dbl x dbl y dbl z\n", unit.c_str()); + for (int i = 0; i < n; i++) + out.Printf(" 0.0 0.0 0.0\n"); + out.CloseFile(); + mprintf("\tWrote '%s' (unit %s, %i atoms).\n", fname.c_str(), unit.c_str(), n); + return 0; +} + +/** Write the remapped structure in the requested format (lib / mol2 / pdb). */ +int Exec_TIMap::WriteAlignedOut(std::string const& fname, AlignedOutFmt fmt, + Topology& top, Frame const& frm, DataSetList const& dsl) const +{ + std::string unit = OffUnitName(ResNameOf(top, "TGT").Truncated(), 'T'); + top.SetParmName(unit, FileName(fname)); + switch (fmt) { + case AOUT_MOL2: + return WriteMol2File(fname, top, frm, dsl); + case AOUT_PDB: + return WritePdbFile(fname, top, frm, dsl); + case AOUT_LIB: + default: + return WriteAmberLib(fname, unit, top, frm); + } +} + +int Exec_TIMap::WriteScmask(std::string const& fname, + Topology const& top0, Topology const& top1, + std::vector const& dual) const +{ + CpptrajFile out; + if (out.OpenWrite(fname)) { + mprinterr("Error: Could not open '%s'\n", fname.c_str()); + return 1; + } + int nShared = 0, nD0 = 0, nD1 = 0; + for (size_t i = 0; i < dual.size(); i++) { + if (dual[i].IsShared()) nShared++; + else if (dual[i].IsTgtOnly()) nD0++; + else nD1++; + } + out.Printf("# timap dual-topology TI masks\n"); + out.Printf("# n_dual= %zu n_shared= %i dummy_in_lambda0= %i dummy_in_lambda1= %i\n", + dual.size(), nShared, nD0, nD1); + out.Printf("# lambda 0 = template (real unmatched, dummy insertions)\n"); + out.Printf("# lambda 1 = target (real insertions, dummy unmatched)\n"); + out.Printf("# pmemd: scmask1 is unique to lambda 0; scmask2 is unique to lambda 1.\n"); + std::string m1, m2; + for (int i = 0; i < (int)dual.size(); i++) { + if (dual[i].IsTplOnly()) { + if (!m1.empty()) m1 += ","; + m1 += ":1@"; + m1 += top0[i].Name().Truncated(); + } + if (dual[i].IsTgtOnly()) { + if (!m2.empty()) m2 += ","; + m2 += ":1@"; + m2 += top1[i].Name().Truncated(); + } + } + out.Printf("scmask1=\"%s\"\n", m1.c_str()); + out.Printf("scmask2=\"%s\"\n", m2.c_str()); + out.CloseFile(); + mprintf("\tWrote '%s'\n", fname.c_str()); + mprintf("\tscmask1 (lambda 0 unique) = %s\n", m1.empty() ? "(none)" : m1.c_str()); + mprintf("\tscmask2 (lambda 1 unique) = %s\n", m2.empty() ? "(none)" : m2.c_str()); + return 0; +} + +/** Per-slot dual-topology table: kind, names, and charges in both end states. */ +int Exec_TIMap::WriteDualAtoms(std::string const& fname, + Topology const& top0, Topology const& top1, + std::vector const& dual) const +{ + CpptrajFile out; + if (out.OpenWrite(fname)) { + mprinterr("Error: Could not open '%s'\n", fname.c_str()); + return 1; + } + out.Printf("# timap dual-topology atoms\n"); + out.Printf("# Kind: SHARED = real in both; TPL_ONLY = dummy in lambda 1; " + "TGT_ONLY = dummy in lambda 0\n"); + out.Printf("%-4s %-8s %-8s %10s %-8s %10s\n", + "#Idx", "Kind", "Name0", "Q0", "Name1", "Q1"); + for (int i = 0; i < (int)dual.size(); i++) { + const char* kind = "SHARED"; + if (dual[i].IsTplOnly()) kind = "TPL_ONLY"; + else if (dual[i].IsTgtOnly()) kind = "TGT_ONLY"; + out.Printf(" %3i %-8s %-8s %10.6f %-8s %10.6f\n", + i + 1, kind, + top0[i].Name().Truncated().c_str(), top0[i].Charge(), + top1[i].Name().Truncated().c_str(), top1[i].Charge()); + } + out.CloseFile(); + mprintf("\tWrote '%s'\n", fname.c_str()); + return 0; +} + +/** Write the human-readable correspondence used by Test_TIMap. */ +int Exec_TIMap::WriteMapFile(std::string const& fname, Topology const& tgt, Topology const& tpl, + TIMatch::Result const& R, std::string const& tgtName, + std::string const& tplName) const +{ + CpptrajFile out; + if (out.OpenWrite(fname)) { + mprinterr("Error: Could not open map file '%s'\n", fname.c_str()); + return 1; + } + out.Printf("# timap tgt='%s' template='%s'\n", tgtName.c_str(), tplName.c_str()); + out.Printf("# kind tgt=%s template=%s\n", R.tgtKind_.c_str(), R.tplKind_.c_str()); + out.Printf("# mapped= %i insertion= %i unmapped_template= %i n_tgt= %i n_tpl= %i\n", + R.nMapped_, R.nInsertion_, R.nUnmappedTpl_, tgt.Natom(), tpl.Natom()); + int aWidth = std::max(6, DigitWidth(tgt.Natom())); + aWidth = std::max(aWidth, DigitWidth(tpl.Natom())); + int nWidth = 6; + for (int i = 0; i < tgt.Natom(); i++) + nWidth = std::max(nWidth, tgt[i].Name().len()); + for (int i = 0; i < tpl.Natom(); i++) + nWidth = std::max(nWidth, tpl[i].Name().len()); + out.Printf("%-*s %*s %-*s %*s %-*s %s\n", + aWidth, "#Out", + aWidth, "TgtAt", + nWidth, "TgtName", + aWidth, "TplAt", + nWidth, "TplName", + "Ins"); + for (int k = 0; k < (int)R.outputOrder_.size(); k++) { + int oldat = R.outputOrder_[k]; + int tplidx = R.mapping_[oldat]; + if (tplidx < 0) + out.Printf("%*i %*i %-*s %*s %-*s %s\n", + aWidth, k + 1, + aWidth, oldat + 1, + nWidth, tgt[oldat].c_str(), + aWidth, "---", + nWidth, "---", + "1"); + else + out.Printf("%*i %*i %-*s %*i %-*s %s\n", + aWidth, k + 1, + aWidth, oldat + 1, + nWidth, tgt[oldat].c_str(), + aWidth, tplidx + 1, + nWidth, tpl[tplidx].c_str(), + "0"); + } + if (R.nUnmappedTpl_ > 0) { + out.Printf("# unmapped template atoms:\n"); + std::vector used(tpl.Natom(), 0); + for (int n = 0; n < tgt.Natom(); n++) { + if (R.mapping_[n] >= 0) used[R.mapping_[n]] = 1; + } + for (int r = 0; r < tpl.Natom(); r++) { + if (!used[r]) + out.Printf("# %i %s\n", r + 1, tpl[r].c_str()); + } + } + out.CloseFile(); + return 0; +} + +/** Filename-safe token from a dataset name (FLE[FLE] -> FLE_FLE). */ +std::string Exec_TIMap::SafeFileToken(std::string const& in) const +{ + std::string out; + for (size_t i = 0; i < in.size(); i++) { + char c = in[i]; + if (std::isalnum((unsigned char)c) || c == '_' || c == '-' || c == '.') + out += c; + else + out += '_'; + } + if (out.empty()) out = "mol"; + return out; +} + +/** Append comma-split tokens from one user string into names. */ +void Exec_TIMap::PushSeriesName(std::vector& names, std::string const& raw) const +{ + if (raw.empty()) return; + if (raw.find(',') == std::string::npos) { + names.push_back(raw); + return; + } + ArgList parts(raw, ","); + for (int i = 0; i < parts.Nargs(); i++) { + std::string s = parts[i]; + // Trim spaces + size_t b = 0; + while (b < s.size() && (s[b] == ' ' || s[b] == '\t')) b++; + size_t e = s.size(); + while (e > b && (s[e - 1] == ' ' || s[e - 1] == '\t')) e--; + if (e > b) names.push_back(s.substr(b, e - b)); + } +} + +/** + * leadopt-style pairwise similarity in [0,1]. + * score = exp(-BETA * delta), BETA = 0.1, with delta = insertions + + * unmatched template atoms (Redmann & Summa leadopt similarity; the + * partial map stands in for a maximum common substructure). + */ +double Exec_TIMap::LeadoptSimScore(int nInsertion, int nUnmappedTpl) const +{ + const double BETA = 0.1; + return std::exp(-BETA * (double)(nInsertion + nUnmappedTpl)); +} + +/** + * Choose the series parent: maximize total mapped atoms when every other + * member is mapped onto the candidate ("shares the most"); on ties prefer + * the minimum structure (fewest atoms), then the similarity sum. + * \return parent index, or -1 on failure. + */ +int Exec_TIMap::SelectSeriesParent(TIMatch const& matcher, + std::vector const& mem, + std::vector& outSimSum, + std::vector& outMappedSum) const +{ + int n = (int)mem.size(); + outSimSum.assign(n, 0.0); + outMappedSum.assign(n, 0); + if (n < 2) return -1; + int best = -1; + for (int i = 0; i < n; i++) { + double simSum = 0.0; + int mappedSum = 0; + for (int j = 0; j < n; j++) { + if (j == i) continue; + TIMatch::Result R; + if (matcher.Match(*mem[j].top, *mem[i].top, R)) { + mprinterr("Error: timap: series match of '%s' onto candidate '%s' failed.\n", + mem[j].name.c_str(), mem[i].name.c_str()); + return -1; + } + mappedSum += R.nMapped_; + simSum += LeadoptSimScore(R.nInsertion_, R.nUnmappedTpl_); + } + outSimSum[i] = simSum; + outMappedSum[i] = mappedSum; + if (best < 0) { + best = i; + continue; + } + bool better = false; + if (mappedSum > outMappedSum[best]) + better = true; + else if (mappedSum == outMappedSum[best]) { + if (mem[i].top->Natom() < mem[best].top->Natom()) + better = true; + else if (mem[i].top->Natom() == mem[best].top->Natom() && + simSum > outSimSum[best] + 1.0e-12) + better = true; + } + if (better) best = i; + } + return best; +} + +/** Parse arguments, run TIMatch::Match, write an aligned structure (lib/mol2/pdb), + * and optionally dual-topology TI files. series mode auto-selects a parent. + */ +// Exec_TIMap::Execute() +Exec::RetType Exec_TIMap::Execute(CpptrajState& State, ArgList& argIn) { + std::string mapout = argIn.GetStringKey("mapout"); + std::string libout = argIn.GetStringKey("out"); + std::string newname = argIn.GetStringKey("name"); + std::string tplName = argIn.GetStringKey("template"); + std::string seedStr = argIn.GetStringKey("seed"); + std::string anchor = argIn.GetStringKey("anchor"); + std::string tiout = argIn.GetStringKey("tiout"); + std::string outprefix = argIn.GetStringKey("outprefix"); + std::string tioutprefix = argIn.GetStringKey("tioutprefix"); + std::string mapoutprefix = argIn.GetStringKey("mapoutprefix"); + std::string parentout = argIn.GetStringKey("parentout"); + std::string fmtStr = argIn.GetStringKey("fmt"); + bool maponly = argIn.hasKey("maponly"); + bool replace = argIn.hasKey("replace"); + bool naorder = argIn.hasKey("naorder"); + bool aaorder = argIn.hasKey("aaorder"); + bool series = argIn.hasKey("series"); + if (naorder && aaorder) { + mprinterr("Error: timap: specify naorder or aaorder, not both.\n"); + return CpptrajState::ERR; + } + AlignedOutFmt outFmt = AOUT_LIB; + if (!fmtStr.empty() && ParseAlignedOutFmt(fmtStr, outFmt)) { + mprinterr("Error: timap: unrecognized fmt '%s' (lib|mol2|pdb).\n", fmtStr.c_str()); + return CpptrajState::ERR; + } + if (!libout.empty()) + outFmt = FmtFromFilename(libout, outFmt); + + // Collect remaining unmarked names (series members, or tgt [/ template]). + std::vector names; + std::string next = argIn.GetStringNext(); + while (!next.empty()) { + PushSeriesName(names, next); + next = argIn.GetStringNext(); + } + // Legacy: template already consumed; if no series and template empty, + // second positional was historically the optional template without the keyword. + // That form is: timap — already collected into names. + + if (series) { + if (!tplName.empty()) { + mprinterr("Error: timap: series chooses the parent; do not also give template.\n"); + return CpptrajState::ERR; + } + if (replace || !newname.empty()) { + mprinterr("Error: timap: series does not support replace / name.\n"); + return CpptrajState::ERR; + } + if (names.size() < 2) { + mprinterr("Error: timap: series needs at least two topologies.\n"); + return CpptrajState::ERR; + } + + TIMatch matcher; + matcher.SetDebug(State.Debug()); + matcher.SetUseNaOrder(naorder); + matcher.SetUseAaOrder(aaorder); + if (aaorder && anchor.empty()) + matcher.SetAnchorName("C"); + if (!anchor.empty()) matcher.SetAnchorName(anchor); + if (!seedStr.empty()) { + std::string l = ToLower(seedStr); + if (l != "auto" && l != "names" && l != "na" && l != "aa" && l != "none") { + mprinterr("Error: timap: unrecognized seed '%s'\n", seedStr.c_str()); + return CpptrajState::ERR; + } + matcher.SetSeed(TIMatch::SeedFromString(seedStr)); + } + + std::vector mem; + mem.reserve(names.size()); + for (size_t i = 0; i < names.size(); i++) { + SeriesMember m; + m.name = names[i]; + m.ds = 0; + m.top = ResolveTop(State, m.name, &m.ds); + if (m.top == 0) { + mprinterr("Error: timap: series member '%s' not found.\n", m.name.c_str()); + return CpptrajState::ERR; + } + mem.push_back(m); + } + + mprintf(" TIMAP: Series of %zu topologies; selecting parent.\n", mem.size()); + mprintf("\tSeed: %s\n", TIMatch::SeedStr( + seedStr.empty() ? TIMatch::SEED_AUTO + : TIMatch::SeedFromString(seedStr))); + if (naorder) + mprintf("\tOrder: nucleic-acid canonical walk of the selected parent.\n"); + else if (aaorder) + mprintf("\tOrder: amino-acid canonical walk (ff19SB) of the selected parent.\n"); + else + mprintf("\tOrder: selected parent's current atom order.\n"); + mprintf("\tParent score: maximize total mapped atoms; ties prefer fewer atoms,\n" + "\t then leadopt-style exp(-0.1*(ins+unmap)) (see util/leadopt citation).\n"); + + std::vector simSum; + std::vector mappedSum; + int pIdx = SelectSeriesParent(matcher, mem, simSum, mappedSum); + if (pIdx < 0) return CpptrajState::ERR; + + SeriesMember const& parent = mem[pIdx]; + mprintf("\tSelected parent: '%s' (%i atoms, sim-sum=%.6f, mapped-sum=%i).\n", + parent.name.c_str(), parent.top->Natom(), simSum[pIdx], mappedSum[pIdx]); + for (size_t i = 0; i < mem.size(); i++) { + mprintf("\t candidate %-16s atoms=%4i sim-sum=%.6f mapped-sum=%i%s\n", + mem[i].name.c_str(), mem[i].top->Natom(), + simSum[i], mappedSum[i], + ((int)i == pIdx) ? " <-- parent" : ""); + } + + if (!parentout.empty()) { + CpptrajFile pout; + if (pout.OpenWrite(parentout)) { + mprinterr("Error: timap: could not write parentout '%s'\n", parentout.c_str()); + return CpptrajState::ERR; + } + pout.Printf("# timap series parent selection\n"); + pout.Printf("# score = sum_j exp(-0.1*(insertions+unmapped)) mapping j onto candidate\n"); + pout.Printf("# adapted from util/leadopt/similarity.py (Redmann & Summa / Summa Lab)\n"); + pout.Printf("parent %s\n", parent.name.c_str()); + pout.Printf("#Idx Name Natom SimSum MappedSum\n"); + for (size_t i = 0; i < mem.size(); i++) { + pout.Printf(" %3zu %-20s %5i %10.6f %10i%s\n", + i + 1, mem[i].name.c_str(), mem[i].top->Natom(), + simSum[i], mappedSum[i], + ((int)i == pIdx) ? " parent" : ""); + } + pout.CloseFile(); + mprintf("\tWrote parent selection to '%s'\n", parentout.c_str()); + } + + bool writeLib = !maponly; + for (size_t i = 0; i < mem.size(); i++) { + SeriesMember const& tgt = mem[i]; + TIMatch::Result R; + if (matcher.Match(*tgt.top, *parent.top, R)) return CpptrajState::ERR; + + mprintf("\t[%zu/%zu] '%s' onto parent '%s': mapped %i / %i " + "(%i insertions, %i parent unmatched).\n", + i + 1, mem.size(), tgt.name.c_str(), parent.name.c_str(), + R.nMapped_, tgt.top->Natom(), R.nInsertion_, R.nUnmappedTpl_); + + std::string tok = SafeFileToken(tgt.name); + std::string thisMap = mapoutprefix.empty() ? std::string() + : mapoutprefix + tok + ".map"; + if (!mapoutprefix.empty()) { + if (WriteMapFile(thisMap, *tgt.top, *parent.top, R, tgt.name, parent.name)) + return CpptrajState::ERR; + mprintf("\t Map written to '%s'\n", thisMap.c_str()); + } + + // TI dual files for every non-parent member (or all if tioutprefix set). + if (!tioutprefix.empty() && (int)i != pIdx) { + std::string tip = tioutprefix + tok; + Frame tgtX, tplX; + LoadCoords(*tgt.top, tgt.ds, tgtX); + LoadCoords(*parent.top, parent.ds, tplX); + Topology top0, top1; + Frame frm0, frm1; + NameType rn0 = ResNameOf(*parent.top, "L0"); + NameType rn1 = ResNameOf(*tgt.top, "L1"); + if (BuildTiUnit(true, *tgt.top, *parent.top, tgtX, tplX, R.dual_, top0, frm0, rn0) || + BuildTiUnit(false, *tgt.top, *parent.top, tgtX, tplX, R.dual_, top1, frm1, rn1)) + { + mprinterr("Error: timap: failed to build TI units for '%s'.\n", tgt.name.c_str()); + return CpptrajState::ERR; + } + std::string f0m = tip + ".0.mol2"; + std::string f1m = tip + ".1.mol2"; + std::string f0l = tip + ".0.lib"; + std::string f1l = tip + ".1.lib"; + std::string fsc = tip + ".scmask"; + std::string fat = tip + ".atoms"; + std::string u0 = OffUnitName(rn0.Truncated(), '0'); + std::string u1 = OffUnitName(rn1.Truncated(), '1'); + if (u0 == u1) { u0 += "0"; u1 += "1"; } + top0.SetParmName(u0, FileName(f0m)); + top1.SetParmName(u1, FileName(f1m)); + if (WriteMol2File(f0m, top0, frm0, State.DSL())) return CpptrajState::ERR; + if (WriteMol2File(f1m, top1, frm1, State.DSL())) return CpptrajState::ERR; + if (WriteAmberLib(f0l, u0, top0, frm0)) return CpptrajState::ERR; + if (WriteAmberLib(f1l, u1, top1, frm1)) return CpptrajState::ERR; + if (WriteScmask(fsc, top0, top1, R.dual_)) return CpptrajState::ERR; + if (WriteDualAtoms(fat, top0, top1, R.dual_)) return CpptrajState::ERR; + if (mapoutprefix.empty()) { + std::string fmap = tip + ".map"; + if (WriteMapFile(fmap, *tgt.top, *parent.top, R, tgt.name, parent.name)) + return CpptrajState::ERR; + } + mprintf("\t TI dual-topology prefix: %s\n", tip.c_str()); + } + + if (!writeLib) continue; + Topology* aligned = tgt.top->ModifyByMap(R.outputOrder_); + if (aligned == 0) { + mprinterr("Error: timap: failed to apply atom order for '%s'.\n", tgt.name.c_str()); + return CpptrajState::ERR; + } + std::string thisOut = outprefix + tok + AlignedOutExt(outFmt); + Frame tgtX, alignedX; + LoadCoords(*tgt.top, tgt.ds, tgtX); + ApplyOrderToFrame(tgtX, R.outputOrder_, alignedX); + if (WriteAlignedOut(thisOut, outFmt, *aligned, alignedX, State.DSL())) { + delete aligned; + return CpptrajState::ERR; + } + mprintf("\t Aligned %s: %s\n", AlignedOutLabel(outFmt), thisOut.c_str()); + delete aligned; + } + return CpptrajState::OK; + } + + // ----- Single-pair / self-map path ----- + if (names.empty()) { + mprinterr("Error: timap: no target topology specified.\n"); + return CpptrajState::ERR; + } + std::string tgtName = names[0]; + if (tplName.empty() && names.size() > 1) + tplName = names[1]; + if (names.size() > 2) { + mprinterr("Error: timap: extra arguments; use 'series' for more than one analog.\n"); + return CpptrajState::ERR; + } + + DataSet* tgtDs = 0; + Topology* tgt = ResolveTop(State, tgtName, &tgtDs); + if (tgt == 0) { + mprinterr("Error: timap: target '%s' not found.\n", tgtName.c_str()); + return CpptrajState::ERR; + } + Topology* tpl = tgt; + DataSet* tplDs = tgtDs; + std::string tplUsed = tgtName; + if (!tplName.empty()) { + tpl = ResolveTop(State, tplName, &tplDs); + if (tpl == 0) { + mprinterr("Error: timap: template '%s' not found.\n", tplName.c_str()); + return CpptrajState::ERR; + } + tplUsed = tplName; + } else { + tplName = tgtName; + tplUsed = tgtName; + } + + TIMatch matcher; + matcher.SetDebug(State.Debug()); + matcher.SetUseNaOrder(naorder); + matcher.SetUseAaOrder(aaorder); + if (aaorder && anchor.empty()) + matcher.SetAnchorName("C"); + if (!anchor.empty()) matcher.SetAnchorName(anchor); + if (!seedStr.empty()) { + std::string l = ToLower(seedStr); + if (l != "auto" && l != "names" && l != "na" && l != "aa" && l != "none") { + mprinterr("Error: timap: unrecognized seed '%s'\n", seedStr.c_str()); + return CpptrajState::ERR; + } + matcher.SetSeed(TIMatch::SeedFromString(seedStr)); + } + + bool writeLib = !libout.empty() || !maponly; + if (writeLib && libout.empty()) { + libout = OffUnitName(ResNameOf(*tgt, "TGT").Truncated(), 'T'); + libout += AlignedOutExt(outFmt); + } + + const char* walkMsg = ""; + if (naorder) walkMsg = " with nucleic-acid canonical walk"; + else if (aaorder) walkMsg = " with amino-acid canonical walk"; + bool selfMap = (tgt == tpl); + if (selfMap) + mprintf(" TIMAP: Reordering '%s' (%i atoms)%s.\n", + tgtName.c_str(), tgt->Natom(), walkMsg); + else + mprintf(" TIMAP: Aligning '%s' (%i atoms) to template '%s' (%i atoms).\n", + tgtName.c_str(), tgt->Natom(), tplUsed.c_str(), tpl->Natom()); + mprintf("\tSeed: %s\n", TIMatch::SeedStr( + seedStr.empty() ? TIMatch::SEED_AUTO + : TIMatch::SeedFromString(seedStr))); + if (naorder) + mprintf("\tOrder: nucleic-acid canonical walk%s.\n", + selfMap ? " of this residue" : " of the template"); + else if (aaorder) + mprintf("\tOrder: amino-acid canonical walk (ff19SB)%s.\n", + selfMap ? " of this residue" : " of the template"); + else if (selfMap) + mprintf("\tOrder: this residue's current atom order (no template, no walk).\n"); + else + mprintf("\tOrder: template file atom order.\n"); + mprintf("\tLeftover insertion anchor: %s\n", + !anchor.empty() ? anchor.c_str() : (aaorder ? "C" : "O3'")); + if (!tiout.empty()) + mprintf("\tTI dual-topology prefix: %s\n", tiout.c_str()); + if (writeLib) + mprintf("\tAligned %s: %s\n", AlignedOutLabel(outFmt), libout.c_str()); + else if (maponly) + mprintf("\tmaponly: not writing an aligned structure file.\n"); + + TIMatch::Result R; + if (matcher.Match(*tgt, *tpl, R)) return CpptrajState::ERR; + + mprintf("\tMapped %i / %i target atoms (%i insertions, %i template atoms unmatched).\n", + R.nMapped_, tgt->Natom(), R.nInsertion_, R.nUnmappedTpl_); + mprintf("\tTarget classified as '%s'; template classified as '%s'.\n", + R.tgtKind_.c_str(), R.tplKind_.c_str()); + if (tgt->Natom() > 0 && R.outputOrder_.size() > 0) { + int first = R.outputOrder_.front(); + int last = R.outputOrder_.back(); + mprintf("\tAligned order starts at %s, ends at %s.\n", + (*tgt)[first].c_str(), (*tgt)[last].c_str()); + } + if (R.nMapped_ == 0) + mprintf("Warning: No atoms were mapped. Check that template and target are related.\n"); + + if (!mapout.empty()) { + mprintf("\tMap written to '%s'\n", mapout.c_str()); + if (WriteMapFile(mapout, *tgt, *tpl, R, tgtName, tplUsed)) + return CpptrajState::ERR; + } + + if (!tiout.empty()) { + if (R.dual_.empty()) { + mprinterr("Error: timap: empty dual-topology layout.\n"); + return CpptrajState::ERR; + } + Frame tgtX, tplX; + bool haveTgtX = LoadCoords(*tgt, tgtDs, tgtX); + bool haveTplX = LoadCoords(*tpl, tplDs, tplX); + if (!haveTgtX || !haveTplX) + mprintf("Warning: Missing coordinates for %s%s%s; dummy/real positions may be zero.\n", + haveTgtX ? "" : "target", + (!haveTgtX && !haveTplX) ? " and " : "", + haveTplX ? "" : "template"); + + Topology top0, top1; + Frame frm0, frm1; + NameType rn0 = ResNameOf(*tpl, "L0"); + NameType rn1 = ResNameOf(*tgt, "L1"); + if (BuildTiUnit(true, *tgt, *tpl, tgtX, tplX, R.dual_, top0, frm0, rn0) || + BuildTiUnit(false, *tgt, *tpl, tgtX, tplX, R.dual_, top1, frm1, rn1)) + { + mprinterr("Error: timap: failed to build TI units.\n"); + return CpptrajState::ERR; + } + if (top0.Natom() != top1.Natom()) { + mprinterr("Error: timap: lambda-0/1 atom counts differ (%i vs %i).\n", + top0.Natom(), top1.Natom()); + return CpptrajState::ERR; + } + mprintf("\tDual topology: %i atoms (%i shared, %i dummy in lambda 0, %i dummy in lambda 1).\n", + top0.Natom(), R.nMapped_, R.nInsertion_, R.nUnmappedTpl_); + + std::string f0m = tiout + ".0.mol2"; + std::string f1m = tiout + ".1.mol2"; + std::string f0l = tiout + ".0.lib"; + std::string f1l = tiout + ".1.lib"; + std::string fsc = tiout + ".scmask"; + std::string fat = tiout + ".atoms"; + std::string u0 = OffUnitName(rn0.Truncated(), '0'); + std::string u1 = OffUnitName(rn1.Truncated(), '1'); + if (u0 == u1) { + u0 += "0"; + u1 += "1"; + } + top0.SetParmName(u0, FileName(f0m)); + top1.SetParmName(u1, FileName(f1m)); + if (WriteMol2File(f0m, top0, frm0, State.DSL())) return CpptrajState::ERR; + if (WriteMol2File(f1m, top1, frm1, State.DSL())) return CpptrajState::ERR; + if (WriteAmberLib(f0l, u0, top0, frm0)) return CpptrajState::ERR; + if (WriteAmberLib(f1l, u1, top1, frm1)) return CpptrajState::ERR; + if (WriteScmask(fsc, top0, top1, R.dual_)) return CpptrajState::ERR; + if (WriteDualAtoms(fat, top0, top1, R.dual_)) return CpptrajState::ERR; + if (mapout.empty()) { + std::string fmap = tiout + ".map"; + mprintf("\tMap written to '%s'\n", fmap.c_str()); + if (WriteMapFile(fmap, *tgt, *tpl, R, tgtName, tplUsed)) + return CpptrajState::ERR; + } + } + + if (!writeLib && newname.empty() && !replace) + return CpptrajState::OK; + + Topology* aligned = tgt->ModifyByMap(R.outputOrder_); + if (aligned == 0) { + mprinterr("Error: timap: failed to apply atom order.\n"); + return CpptrajState::ERR; + } + if (writeLib) { + Frame tgtX, alignedX; + LoadCoords(*tgt, tgtDs, tgtX); + ApplyOrderToFrame(tgtX, R.outputOrder_, alignedX); + if (WriteAlignedOut(libout, outFmt, *aligned, alignedX, State.DSL())) { + delete aligned; + return CpptrajState::ERR; + } + } + + if (replace) { + if (tgtDs != 0 && tgtDs->Type() == DataSet::TOPOLOGY) { + ((DataSet_Topology*)tgtDs)->SetTop(*aligned); + mprintf("\tReplaced topology '%s' with aligned atom order.\n", tgtName.c_str()); + } else { + mprinterr("Error: timap: replace requires to be a topology set.\n"); + delete aligned; + return CpptrajState::ERR; + } + } + if (!newname.empty()) { + aligned->SetParmName(newname, tgt->OriginalFilename()); + if (State.AddTopology(*aligned, newname)) { + delete aligned; + return CpptrajState::ERR; + } + mprintf("\tAligned topology added as '%s'\n", newname.c_str()); + } + delete aligned; + return CpptrajState::OK; +} diff --git a/src/Exec_TIMap.h b/src/Exec_TIMap.h new file mode 100644 index 0000000000..b2eecbb22b --- /dev/null +++ b/src/Exec_TIMap.h @@ -0,0 +1,73 @@ +#ifndef INC_EXEC_TIMAP_H +#define INC_EXEC_TIMAP_H +#include "Exec.h" +#include "NameType.h" +#include "TIMatch.h" +#include +#include +class Atom; +class DataSet; +class DataSetList; +class FileName; +class Frame; +class NameType; +class Topology; +/// Immediate command: write a residue library in template atom order; optional TI pair. +/** Command: timap (alias: timatch). See Help() and 'help timap extended'. + * \author Nathan D. Levinzon + */ +class Exec_TIMap : public Exec { + public: + Exec_TIMap() : Exec(PARM) {} + void Help() const; + void Help(ArgList&) const; + DispatchObject* Alloc() const { return (DispatchObject*)new Exec_TIMap(); } + RetType Execute(CpptrajState&, ArgList&); + private: + /// Aligned-structure output format (default Amber OFF .lib for LEaP/TI). + enum AlignedOutFmt { AOUT_LIB = 0, AOUT_MOL2, AOUT_PDB }; + /// One topology in a series-mode parent search. + struct SeriesMember { + std::string name; + Topology* top; + DataSet* ds; + }; + + Topology* FindNamedTop(DataSetList&, std::string const&, DataSet**) const; + void ParseTopSpec(std::string const&, std::string&, std::string&) const; + std::string FileStem(FileName const&) const; + Topology* ResolveTop(CpptrajState&, std::string const&, DataSet**) const; + bool LoadCoords(Topology const&, DataSet*, Frame&) const; + void CopyXyz(Frame const&, int, Frame&, int) const; + void ApplyOrderToFrame(Frame const&, std::vector const&, Frame&) const; + std::string OffUnitName(std::string, char) const; + NameType ResNameOf(Topology const&, const char*) const; + int BuildTiUnit(bool, Topology const&, Topology const&, Frame const&, Frame const&, + std::vector const&, + Topology&, Frame&, NameType const&) const; + int WriteMol2File(std::string const&, Topology&, Frame const&, DataSetList const&) const; + int WritePdbFile(std::string const&, Topology&, Frame const&, DataSetList const&) const; + const char* AlignedOutExt(AlignedOutFmt) const; + const char* AlignedOutLabel(AlignedOutFmt) const; + int ParseAlignedOutFmt(std::string const&, AlignedOutFmt&) const; + AlignedOutFmt FmtFromFilename(std::string const&, AlignedOutFmt) const; + std::string Q(std::string const&) const; + std::string AtomTypeStr(Atom const&) const; + int AtomNum1(Topology const&, const char*) const; + void LeapConnect(Topology const&, int&, int&, const char*&) const; + int WriteAmberLib(std::string const&, std::string const&, Topology const&, Frame const&) const; + int WriteAlignedOut(std::string const&, AlignedOutFmt, Topology&, Frame const&, + DataSetList const&) const; + int WriteScmask(std::string const&, Topology const&, Topology const&, + std::vector const&) const; + int WriteDualAtoms(std::string const&, Topology const&, Topology const&, + std::vector const&) const; + int WriteMapFile(std::string const&, Topology const&, Topology const&, + TIMatch::Result const&, std::string const&, std::string const&) const; + std::string SafeFileToken(std::string const&) const; + void PushSeriesName(std::vector&, std::string const&) const; + double LeadoptSimScore(int, int) const; + int SelectSeriesParent(TIMatch const&, std::vector const&, + std::vector&, std::vector&) const; +}; +#endif diff --git a/src/TIMatch.cpp b/src/TIMatch.cpp new file mode 100644 index 0000000000..89cbd478fc --- /dev/null +++ b/src/TIMatch.cpp @@ -0,0 +1,1681 @@ +#include "TIMatch.h" +#include "Topology.h" +#include "CpptrajStdio.h" +#include "StringRoutines.h" +#include +#include +#include +#include +#include + +/* Source encoding: UTF-8. Unicode (O5′, χ, λ, 2′) is comments-only; string + * literals stay ASCII so this file compiles without a UTF-8 compiler flag. + * + * ============================================================================= + * TIMatch — partial graph match for thermodynamic integration (timap) + * ============================================================================= + * + * PURPOSE + * Build a correspondence between a *target* residue and a *template* + * (parent) residue so shared atoms share indices for TI / dual topology. + * Unlike atommap, a 1:1 map is *not* required: extras on the target are + * insertions; unmatched template atoms are simply skipped in the target + * permutation (and become dummies in DualOrder). + * + * PIPELINE (Match()) + * 1. Graph Bond adjacency for tgt and tpl (all / heavy / H lists). + * 2. DetectScaffold Classify each residue: nucleotide | sugar | base | + * amino | unknown; fill NA roles and/or peptide aa_ roles. + * 3. MapGraphs Seed → Grow → unique names → Grow → hydrogens. + * Result: mapping_[tgt] = tpl index or −1 (insertion). + * 4. ParentOrder Template walk: file order, naorder, or aaorder. + * 5. OutputOrder Permute tgt atoms to follow that walk + insertions. + * 6. DualOrder Same layout but keep unmatched tpl slots (TI dummies). + * + * CHEMISTRIES (same code paths) + * - Nucleic acids: scaffold roles seed (FLE rA vs ERN dA; 2′-OH inserts). + * - Amino acids: N/CA/C/O/CB seed (Cys vs Sec; Se inserts at CB). + * - Small mols: unique names + grow (benzene vs phenol; OH inserts). + * + * KEY INVARIANTS (do not break without updating Exec_TIMap + tests) + * - mapping_[i] is either −1 or a unique template index (no two tgts share). + * - outputOrder_ is a permutation of 0..Ntgt−1 (ModifyByMap requires this). + * - dual_.size() == nMapped + nInsertion + nUnmappedTpl. + * - Grow signatures use *element* + heavy degree + neighbor elements so + * S≠Se and O≠H at the same site become insertions, not forced matches. + * - Unique-name seed also requires matching AtomicNumber/Element. + * + * WHERE TO CHANGE WHAT + * - New NA parent walk: CanonicalWalk / DetectScaffold roles. + * - New AA walk (force field): CanonicalAaWalk / DetectAmino. + * - Matching too greedy/shy: MakeSig, Grow, MatchUniqueNames. + * - Insertion placement: OutputOrder / DualOrder / EmitTree. + * - Dual-topology layout: DualOrder only (OutputOrder skips tpl-only). + * + * Nested Graph / Scaffold types live in an anonymous namespace so free + * helpers can use them without polluting the TIMatch public API. + */ + +// ----------------------------------------------------------------------------- +TIMatch::TIMatch() : + seed_(SEED_AUTO), + debug_(0), + useNaOrder_(false), + useAaOrder_(false), + anchorName_("O3'") +{} + +const char* TIMatch::SeedStr(SeedType t) { + switch (t) { + case SEED_AUTO: return "auto"; + case SEED_NAMES: return "names"; + case SEED_NA: return "na"; + case SEED_AA: return "aa"; + case SEED_NONE: return "none"; + } + return "auto"; +} + +/// Parse seed {auto|names|na|aa|none}; unknown strings fall back to auto. +TIMatch::SeedType TIMatch::SeedFromString(std::string const& s) { + std::string l = ToLower(s); + if (l == "auto") return SEED_AUTO; + if (l == "names") return SEED_NAMES; + if (l == "na") return SEED_NA; + if (l == "aa") return SEED_AA; + if (l == "none") return SEED_NONE; + return SEED_AUTO; +} + +namespace { + +typedef std::vector Iarray; + +// ============================================================================= +// Graph helpers — element predicates and bond adjacency +// ============================================================================= + +/// True for elemental hydrogen (not extra points / Drude particles). +static inline bool IsH(Atom const& a) { + return (a.Element() == Atom::HYDROGEN); +} + +/// Heavy = not H, extra point, or Drude. Used for ring detection and signatures. +static inline bool IsHeavy(Atom const& a) { + return (!IsH(a) && a.Element() != Atom::EXTRAPT && a.Element() != Atom::DRUDE); +} + +/// Truncated Amber name (no padding). Comparisons and map files use this. +static std::string AtomName(Atom const& a) { + return a.Name().Truncated(); +} + +// ----------------------------------------------------------------------------- +/// Bonded graph over one Topology: full, heavy-only, and hydrogen neighbor lists. +/** Built once per residue so DetectScaffold / Grow / CanonicalWalk can ask + * "who is bonded to C2′?" without walking Atom::bond_iterator each time. + * + * Indices are Topology atom indices (0-based). nbr_[i] lists every bonded + * partner; heavyNbr_ / hNbr_ are filtered views for grow signatures and + * hydrogen zip. Topology is not owned — it must outlive the Graph. + */ +class Graph { + public: + Graph() : top_(0) {} + explicit Graph(Topology const& t) { Setup(t); } + + /// Fill nbr_ / heavyNbr_ / hNbr_ from Topology bonds. Indices match atom #. + void Setup(Topology const& t) { + top_ = &t; + int n = t.Natom(); + nbr_.assign(n, Iarray()); + heavyNbr_.assign(n, Iarray()); + hNbr_.assign(n, Iarray()); + for (int i = 0; i < n; i++) { + Atom const& at = t[i]; + for (Atom::bond_iterator b = at.bondbegin(); b != at.bondend(); ++b) { + nbr_[i].push_back(*b); + if (IsHeavy(t[*b])) + heavyNbr_[i].push_back(*b); + else if (IsH(t[*b])) + hNbr_[i].push_back(*b); + } + } + } + + Topology const& Top() const { return *top_; } + int Natom() const { return top_->Natom(); } + Atom const& operator[](int i) const { return (*top_)[i]; } + std::string Name(int i) const { return AtomName((*top_)[i]); } + Atom::AtomicElementType Elt(int i) const { return (*top_)[i].Element(); } + bool IsHydrogen(int i) const { return IsH((*top_)[i]); } + int HeavyDeg(int i) const { return (int)heavyNbr_[i].size(); } + Iarray const& Nbr(int i) const { return nbr_[i]; } + Iarray const& HeavyNbr(int i) const { return heavyNbr_[i]; } + Iarray const& HNbr(int i) const { return hNbr_[i]; } + + /// Every atom whose truncated name equals \p name (0, 1, or many hits). + Iarray FindByName(std::string const& name) const { + Iarray hits; + for (int i = 0; i < Natom(); i++) { + if (Name(i) == name) + hits.push_back(i); + } + return hits; + } + + Iarray HeavyNbrsOf(int i) const { return heavyNbr_[i]; } + + private: + Topology const* top_; + std::vector nbr_; + std::vector heavyNbr_; + std::vector hNbr_; +}; + +// ============================================================================= +// Scaffold — chemical roles on one residue (NA and/or peptide) +// ============================================================================= + +/// Nucleic-acid (or ligand) roles detected on one residue. +/** role_[P] … role_[CHI] are atom indices or −1. sugarRing_ is the chosen + * furanose (4 C + 1 O/S). base_ is the connected component hanging off χ, + * excluding C1′. kind_ is nucleotide | sugar | base | amino | unknown. + * + * Amino acids fill aa_ (N, CA, C, O, CB) and kind_ "amino" when the + * peptide N–CA–C=O motif is found. Small molecules stay kind "unknown". + * + * ok_ means "looks like a nucleotide scaffold good enough to SEED_NA". + * aaOk_ means backbone N/CA/C were found (enough for SEED_AA). + * notes_ carries soft failures for debug / map headers (never fatal alone). + */ +class Scaffold { + public: + /// Canonical nucleotide roles. CHI is the glycosidic N (or C) of χ. + enum Role { P = 0, O5p, C5p, C4p, O4p, C1p, C2p, C3p, O3p, CHI, NROLES }; + /// Peptide backbone / β-carbon roles (ff19SB naming). + enum AaRole { AA_N = 0, AA_CA, AA_C, AA_O, AA_CB, AA_NROLES }; + static const char* RoleStr(Role r) { + static const char* n[NROLES] = { + "P", "O5'", "C5'", "C4'", "O4'", "C1'", "C2'", "C3'", "O3'", "chi" + }; + return n[r]; + } + static const char* AaRoleStr(AaRole r) { + static const char* n[AA_NROLES] = { "N", "CA", "C", "O", "CB" }; + return n[r]; + } + + Scaffold() : hasP_(false), ok_(false), aaOk_(false), isPro_(false) { + for (int i = 0; i < NROLES; i++) role_[i] = -1; + for (int i = 0; i < AA_NROLES; i++) aa_[i] = -1; + kind_ = "unknown"; + family_ = "none"; + twoPrime_ = "unknown"; + } + + int Get(Role r) const { return role_[r]; } + void Set(Role r, int i) { role_[r] = i; } + int Aa(AaRole r) const { return aa_[r]; } + void SetAa(AaRole r, int i) { aa_[r] = i; } + bool HasSugar() const { return !sugarRing_.empty(); } + /// Combined nucleotide: furanose + glycosidic N + a nucleobase component. + bool IsCombined() const { + return HasSugar() && role_[CHI] >= 0 && !base_.empty(); + } + + int role_[NROLES]; + Iarray op_; ///< Non-bridging phosphate oxygens (name-sorted). + Iarray sugarRing_; ///< Atoms of the chosen furanose (order from cycle search). + Iarray base_; ///< χ-component excluding C1′ (sorted unique indices). + std::string kind_; ///< nucleotide | sugar | base | amino | unknown + std::string family_; ///< purine | pyrimidine | pro | gly | std | none | … + std::string twoPrime_; ///< HH | HOH | OMe | OR | F | other | unknown + bool hasP_; + bool ok_; ///< True if NA scaffold is seed-worthy. + bool aaOk_; ///< True if N, CA, C were assigned. + bool isPro_; ///< N bonded to ≥2 carbons (proline-like). + int aa_[AA_NROLES]; + std::vector notes_; +}; + +// ============================================================================= +// Ring finding — furanose / nucleobase detection +// ============================================================================= + +typedef std::vector Cycles; + +/// Rotate a cycle so the lowest index is first; pick the lexicographically +/// smaller of the two directions so the same ring is not stored twice. +static Iarray NormalizeCycle(Iarray path) { + if (path.empty()) return path; + int i0 = 0; + for (int i = 1; i < (int)path.size(); i++) { + if (path[i] < path[i0]) i0 = i; + } + Iarray rot; + rot.reserve(path.size()); + for (int k = 0; k < (int)path.size(); k++) + rot.push_back(path[(i0 + k) % path.size()]); + if (rot.size() > 2 && rot.back() < rot[1]) { + Iarray rev(1, rot[0]); + for (int k = (int)rot.size() - 1; k >= 1; k--) + rev.push_back(rot[k]); + int j0 = 0; + for (int i = 1; i < (int)rev.size(); i++) { + if (rev[i] < rev[j0]) j0 = i; + } + Iarray rot2; + for (int k = 0; k < (int)rev.size(); k++) + rot2.push_back(rev[(j0 + k) % rev.size()]); + return rot2; + } + return rot; +} + +/// Depth-first search for simple heavy-atom cycles of a fixed length. +/** \p start is the origin of this search (only neighbors ≥ start are + * followed, which together with NormalizeCycle deduplicates rings). + */ +static void DfsCycle(int start, Iarray const& path, std::set const& seen, + std::vector const& heavyNbr, int length, + std::set& found) +{ + if ((int)path.size() == length) { + int last = path.back(); + Iarray const& nbr = heavyNbr[last]; + if (std::find(nbr.begin(), nbr.end(), start) != nbr.end() && path[0] == start) + found.insert(NormalizeCycle(path)); + return; + } + int last = path.back(); + Iarray const& nbr = heavyNbr[last]; + for (Iarray::const_iterator it = nbr.begin(); it != nbr.end(); ++it) { + int nxt = *it; + if (seen.find(nxt) != seen.end()) continue; + if (nxt < start) continue; + Iarray path2 = path; + path2.push_back(nxt); + std::set seen2 = seen; + seen2.insert(nxt); + DfsCycle(start, path2, seen2, heavyNbr, length, found); + } +} + +/// All unique heavy-atom cycles of \p length (5 = furanose / imidazole, 6 = pyrimidine). +static Cycles CyclesOfLength(Graph const& g, int length) { + std::set found; + std::vector heavy; + heavy.reserve(g.Natom()); + for (int i = 0; i < g.Natom(); i++) + heavy.push_back(g.HeavyNbr(i)); + for (int i = 0; i < g.Natom(); i++) { + Iarray path(1, i); + std::set seen; + seen.insert(i); + DfsCycle(i, path, seen, heavy, length, found); + } + Cycles out(found.begin(), found.end()); + return out; +} + +/// Count O / C / N / S among the atoms of one cycle (used to score sugar vs nucleobase rings). +static void CountRingElts(Graph const& g, Iarray const& cyc, + int& nO, int& nC, int& nN, int& nS) +{ + nO = nC = nN = nS = 0; + for (Iarray::const_iterator it = cyc.begin(); it != cyc.end(); ++it) { + switch (g.Elt(*it)) { + case Atom::OXYGEN: nO++; break; + case Atom::CARBON: nC++; break; + case Atom::NITROGEN: nN++; break; + case Atom::SULFUR: nS++; break; + default: break; + } + } +} + +/// Furanose: four carbons and exactly one O or S, no nitrogen (so not imidazole). +static bool IsSugarRing(Graph const& g, Iarray const& cyc) { + int nO, nC, nN, nS; + CountRingElts(g, cyc, nO, nC, nN, nS); + return (nC == 4 && nN == 0 && (nO + nS) == 1); +} + +/// Prefer the 5-ring whose carbons have the most nucleotide-like exo substituents +/// (glycosidic N, C5′ carbon, O3′/O5′ oxygen). Ties break lexicographically. +/** Score weights (per ring carbon): exo N +4, exo C +3, exo O +2. Imidazole + * 5-rings fail IsSugarRing (they contain N) so they never compete here. + */ +static Iarray PickSugarRing(Graph const& g, Cycles const& rings5) { + int bestScore = -1; + Iarray best; + for (Cycles::const_iterator c = rings5.begin(); c != rings5.end(); ++c) { + if (!IsSugarRing(g, *c)) continue; + std::set ring(c->begin(), c->end()); + int score = 0; + for (Iarray::const_iterator it = c->begin(); it != c->end(); ++it) { + if (g.Elt(*it) != Atom::CARBON) continue; + bool exoC = false, exoO = false, exoN = false; + Iarray const& hn = g.HeavyNbr(*it); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (ring.find(*j) != ring.end()) continue; + if (g.Elt(*j) == Atom::CARBON) exoC = true; + if (g.Elt(*j) == Atom::OXYGEN) exoO = true; + if (g.Elt(*j) == Atom::NITROGEN) exoN = true; + } + if (exoC) score += 3; + if (exoO) score += 2; + if (exoN) score += 4; + } + if (score > bestScore || (score == bestScore && (best.empty() || *c < best))) { + bestScore = score; + best = *c; + } + } + return best; +} + +/// Connected component of χ, stopping at C1′ so the sugar is not swallowed into the base. +static Iarray BaseAtoms(Graph const& g, int chi, int c1) { + std::set seen; + seen.insert(chi); + if (c1 >= 0) seen.insert(c1); + Iarray stack(1, chi); + Iarray out; + while (!stack.empty()) { + int i = stack.back(); + stack.pop_back(); + out.push_back(i); + Iarray const& nbr = g.Nbr(i); + for (Iarray::const_iterator j = nbr.begin(); j != nbr.end(); ++j) { + if (seen.find(*j) != seen.end()) continue; + if (c1 >= 0 && *j == c1) continue; + seen.insert(*j); + stack.push_back(*j); + } + } + std::sort(out.begin(), out.end()); + out.erase(std::unique(out.begin(), out.end()), out.end()); + return out; +} + +/// purine if an imidazole 5-ring and a 6-ring with N are both present; else pyrimidine. +static std::string BaseFamily(Graph const& g, int chi, + Cycles const& rings5, Cycles const& rings6) +{ + bool imidazole = false; + for (Cycles::const_iterator c = rings5.begin(); c != rings5.end(); ++c) { + int nO, nC, nN, nS; + CountRingElts(g, *c, nO, nC, nN, nS); + if (nN >= 1 && nO == 0) imidazole = true; + } + bool sixN = false; + for (Cycles::const_iterator c = rings6.begin(); c != rings6.end(); ++c) { + int nO, nC, nN, nS; + CountRingElts(g, *c, nO, nC, nN, nS); + if (nN >= 1) sixN = true; + } + if (imidazole && sixN) return "purine"; + if (sixN) return "pyrimidine"; + if (imidazole) return "purine"; + int nC = 0; + Iarray const& hn = g.HeavyNbr(chi); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (g.Elt(*j) == Atom::CARBON) nC++; + } + if (nC >= 2) return "purine"; + return "pyrimidine"; +} + +/// Classify the 2′ substituent: HH (deoxy), HOH (ribose), OMe, OR, F, other. +static std::string TwoPrimePattern(Graph const& g, + Scaffold const& sc) +{ + int c2 = sc.Get(Scaffold::C2p); + if (c2 < 0) return "unknown"; + std::set skip; + if (sc.Get(Scaffold::C1p) >= 0) skip.insert(sc.Get(Scaffold::C1p)); + if (sc.Get(Scaffold::C3p) >= 0) skip.insert(sc.Get(Scaffold::C3p)); + if (sc.Get(Scaffold::O4p) >= 0) skip.insert(sc.Get(Scaffold::O4p)); + int nO = 0, nC = 0, nF = 0, nH = (int)g.HNbr(c2).size(); + int oIdx = -1; + Iarray const& hn = g.HeavyNbr(c2); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (skip.find(*j) != skip.end()) continue; + if (g.Elt(*j) == Atom::OXYGEN) { nO++; oIdx = *j; } + else if (g.Elt(*j) == Atom::CARBON) nC++; + else if (g.Elt(*j) == Atom::FLUORINE) nF++; + } + if (nO == 0 && nC == 0 && nF == 0 && nH >= 1) return "HH"; + if (nO == 1 && nC == 0 && oIdx >= 0) { + if ((int)g.HeavyNbr(oIdx).size() == 1) return "HOH"; + Iarray carbons; + Iarray const& ohn = g.HeavyNbr(oIdx); + for (Iarray::const_iterator j = ohn.begin(); j != ohn.end(); ++j) { + if (*j != c2 && g.Elt(*j) == Atom::CARBON) + carbons.push_back(*j); + } + if (carbons.size() == 1 && g.HeavyDeg(carbons[0]) == 1) return "OMe"; + if (!carbons.empty()) return "OR"; + return "OR"; + } + if (nF == 1) return "F"; + return "other"; +} + +/// Best glycosidic nitrogen for a nucleobase fragment (no sugar): N9 / N1 names win. +static int GuessGlycosidicN(Graph const& g, + Cycles const& rings5, Cycles const& rings6) +{ + int best = -1; + int bestScore = -999; + for (int i = 0; i < g.Natom(); i++) { + if (g.Elt(i) != Atom::NITROGEN) continue; + int nC = 0; + Iarray const& hn = g.HeavyNbr(i); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (g.Elt(*j) == Atom::CARBON) nC++; + } + bool in5 = false, in6 = false; + for (Cycles::const_iterator c = rings5.begin(); c != rings5.end(); ++c) { + if (std::find(c->begin(), c->end(), i) == c->end()) continue; + int nO, nC2, nN, nS; + CountRingElts(g, *c, nO, nC2, nN, nS); + if (nN >= 1 && nO == 0) in5 = true; + } + for (Cycles::const_iterator c = rings6.begin(); c != rings6.end(); ++c) { + if (std::find(c->begin(), c->end(), i) != c->end()) in6 = true; + } + int score = nC; + if (in5) score += 8; + if (in6 && !in5) score += 3; + std::string nm = ToLower(g.Name(i)); + if (nm == "n9") score += 10; + else if (nm == "n1" && !in5) score += 10; + if (score > bestScore) { + bestScore = score; + best = i; + } + } + return best; +} + +/// Fill any still-empty roles from unique Amber names (P, O5′, C1′, N9, …). +static void NameHints(Graph const& g, Scaffold& sc) { + static const char* names[] = { + "P", "O5'", "C5'", "C4'", "O4'", "C1'", "C2'", "C3'", "O3'", "N9", "N1", "S4'", 0 + }; + static const Scaffold::Role roles[] = { + Scaffold::P, Scaffold::O5p, Scaffold::C5p, + Scaffold::C4p, Scaffold::O4p, Scaffold::C1p, + Scaffold::C2p, Scaffold::C3p, Scaffold::O3p, + Scaffold::CHI, Scaffold::CHI, Scaffold::O4p + }; + for (int k = 0; names[k] != 0; k++) { + Scaffold::Role role = roles[k]; + if (sc.Get(role) >= 0) continue; + Iarray hits = g.FindByName(names[k]); + if (hits.size() == 1) + sc.Set(role, hits[0]); + } +} + +/// If the topology carries ModXNA HEAD/TAIL/ANCHOR metadata, use it as χ / O3′ / C1′. +static void ApplyModxna(Graph const& g, Scaffold& sc) { + ModXNA_Info const& mx = g.Top().Modxna(); + if (!mx.HasModxna()) return; + std::string head = mx.Head(); + if (!head.empty() && head[0] == '@') head = head.substr(1); + std::string tail = mx.Tail(); + if (!tail.empty() && tail[0] == '@') tail = tail.substr(1); + std::string anchor = mx.Anchor(); + if (!anchor.empty() && anchor[0] == '@') anchor = anchor.substr(1); + Iarray hh = head.empty() ? Iarray() : g.FindByName(head); + Iarray tt = tail.empty() ? Iarray() : g.FindByName(tail); + Iarray aa = anchor.empty() ? Iarray() : g.FindByName(anchor); + int hi = (hh.size() == 1) ? hh[0] : -1; + int ti = (tt.size() == 1) ? tt[0] : -1; + int ai = (aa.size() == 1) ? aa[0] : -1; + if (hi >= 0 && sc.Get(Scaffold::CHI) < 0 && + (g.Elt(hi) == Atom::NITROGEN || g.Elt(hi) == Atom::CARBON)) + sc.Set(Scaffold::CHI, hi); + if (ti >= 0 && sc.Get(Scaffold::O3p) < 0 && g.Elt(ti) == Atom::OXYGEN) + sc.Set(Scaffold::O3p, ti); + if (ai >= 0 && sc.Get(Scaffold::C1p) < 0) + sc.Set(Scaffold::C1p, ai); +} + +/// Assign C1′/C4′ (the two carbons on O4′), then C2′/C3′, C5′, O5′, OP, O3′. +/** Scoring for the O4′-bonded carbons (lower score → C1′): + * −10 exo N (glycosidic attachment — strong C1′ signal) + * +6 exo C that itself bears O (C5′/O5′ side — C4′ signal) + * +2 exo O + * −4 exo methyl-like C (degree ≤1) — often a C1′ base-cap carbon + * + * After C1′/C4′: C3′ is the remaining ring C bonded to C4′; C2′ is the other. + * Swap C2′/C3′ if C2′ is not bonded to C1′. C5′ is exo-C on C4′; O5′ prefers + * the oxygen also bonded to P when a phosphate is present. OP = other O on P. + * O3′ = exo-O on C3′. + */ +static void AssignSugarRoles(Graph const& g, Scaffold& sc, + Iarray const& sugar) +{ + std::set ring(sugar.begin(), sugar.end()); + int o4 = -1; + Iarray ringC; + for (Iarray::const_iterator it = sugar.begin(); it != sugar.end(); ++it) { + if (g.Elt(*it) == Atom::OXYGEN || g.Elt(*it) == Atom::SULFUR) o4 = *it; + if (g.Elt(*it) == Atom::CARBON) ringC.push_back(*it); + } + sc.Set(Scaffold::O4p, o4); + Iarray bondedToO; + if (o4 >= 0) { + for (Iarray::const_iterator it = ringC.begin(); it != ringC.end(); ++it) { + Iarray const& hn = g.HeavyNbr(*it); + if (std::find(hn.begin(), hn.end(), o4) != hn.end()) + bondedToO.push_back(*it); + } + } + int c1 = -1, c4 = -1; + if (bondedToO.size() == 2) { + std::vector< std::pair > scored; + for (Iarray::const_iterator it = bondedToO.begin(); it != bondedToO.end(); ++it) { + int score = 0; + Iarray const& hn = g.HeavyNbr(*it); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (ring.find(*j) != ring.end()) continue; + if (g.Elt(*j) == Atom::NITROGEN) score -= 10; + else if (g.Elt(*j) == Atom::OXYGEN) score += 2; + else if (g.Elt(*j) == Atom::CARBON) { + int nO = 0; + Iarray const& jn = g.HeavyNbr(*j); + for (Iarray::const_iterator k = jn.begin(); k != jn.end(); ++k) + if (g.Elt(*k) == Atom::OXYGEN) nO++; + if (nO) score += 6; + else if ((int)jn.size() <= 1) score -= 4; + else score += 1; + } + } + scored.push_back(std::make_pair(score, *it)); + } + std::sort(scored.begin(), scored.end()); + c1 = scored[0].second; + c4 = scored[1].second; + } else if (bondedToO.size() == 1) { + c1 = bondedToO[0]; + } + sc.Set(Scaffold::C1p, c1); + sc.Set(Scaffold::C4p, c4); + if (c1 < 0 || c4 < 0) { + sc.notes_.push_back("could not split C1'/C4' on the sugar ring"); + return; + } + Iarray remain; + for (Iarray::const_iterator it = ringC.begin(); it != ringC.end(); ++it) { + if (*it != c1 && *it != c4) remain.push_back(*it); + } + int c3 = -1; + for (Iarray::const_iterator it = remain.begin(); it != remain.end(); ++it) { + Iarray const& hn = g.HeavyNbr(*it); + if (std::find(hn.begin(), hn.end(), c4) != hn.end()) { c3 = *it; break; } + } + int c2 = -1; + for (Iarray::const_iterator it = remain.begin(); it != remain.end(); ++it) { + if (*it != c3) { c2 = *it; break; } + } + if (c2 >= 0 && c1 >= 0) { + Iarray const& hn = g.HeavyNbr(c2); + bool c2c1 = std::find(hn.begin(), hn.end(), c1) != hn.end(); + if (!c2c1 && c3 >= 0) { + Iarray const& hn3 = g.HeavyNbr(c3); + if (std::find(hn3.begin(), hn3.end(), c1) != hn3.end()) + std::swap(c2, c3); + } + } + sc.Set(Scaffold::C3p, c3); + sc.Set(Scaffold::C2p, c2); + + if (c4 >= 0) { + Iarray exoC; + Iarray const& hn = g.HeavyNbr(c4); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (ring.find(*j) == ring.end() && g.Elt(*j) == Atom::CARBON) + exoC.push_back(*j); + } + if (!exoC.empty()) { + int c5 = exoC[0]; + sc.Set(Scaffold::C5p, c5); + Iarray o5cand; + Iarray const& c5n = g.HeavyNbr(c5); + for (Iarray::const_iterator j = c5n.begin(); j != c5n.end(); ++j) { + if (g.Elt(*j) == Atom::OXYGEN) o5cand.push_back(*j); + } + int p = sc.Get(Scaffold::P); + if (p >= 0) { + int o5 = -1; + for (Iarray::const_iterator j = o5cand.begin(); j != o5cand.end(); ++j) { + Iarray const& jn = g.HeavyNbr(*j); + if (std::find(jn.begin(), jn.end(), p) != jn.end()) { o5 = *j; break; } + } + sc.Set(Scaffold::O5p, o5 >= 0 ? o5 : (o5cand.empty() ? -1 : o5cand[0])); + } else if (!o5cand.empty()) { + sc.Set(Scaffold::O5p, o5cand[0]); + } + } + } + + int p = sc.Get(Scaffold::P); + int o5 = sc.Get(Scaffold::O5p); + if (p >= 0) { + Iarray ops; + Iarray const& pn = g.HeavyNbr(p); + for (Iarray::const_iterator j = pn.begin(); j != pn.end(); ++j) { + if (g.Elt(*j) == Atom::OXYGEN && *j != o5) ops.push_back(*j); + } + std::sort(ops.begin(), ops.end(), [&g](int a, int b){ return g.Name(a) < g.Name(b); }); + sc.op_ = ops; + } + + if (c3 >= 0) { + Iarray o3cand; + Iarray const& hn = g.HeavyNbr(c3); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (ring.find(*j) == ring.end() && g.Elt(*j) == Atom::OXYGEN) + o3cand.push_back(*j); + } + if (!o3cand.empty()) + sc.Set(Scaffold::O3p, o3cand[0]); + } +} + +/// Unique truncated name, or −1 if missing / not unique. +static int UniqueName(Graph const& g, const char* nm) { + Iarray hits = g.FindByName(nm); + return (hits.size() == 1) ? hits[0] : -1; +} + +/// Terminal oxygen (heavy degree 1): carbonyl O, carboxylate O, hydroxyl O. +static bool IsTermO(Graph const& g, int i) { + return (g.Elt(i) == Atom::OXYGEN && g.HeavyDeg(i) == 1); +} + +/// Fill peptide N / CA / C / O / CB from Amber names, then the N–CA–C=O motif. +/** ff19SB amino19.lib is the baseline: N, H, CA, HA, side chain from CB, C, O + * (proline: N, CD … CB, CA, C, O). Noncanonical residues that keep backbone + * names are labeled the same way; graph backup handles renamed backbones. + */ +static void DetectAmino(Graph const& g, Scaffold& sc) { + sc.SetAa(Scaffold::AA_N, UniqueName(g, "N")); + sc.SetAa(Scaffold::AA_CA, UniqueName(g, "CA")); + sc.SetAa(Scaffold::AA_C, UniqueName(g, "C")); + sc.SetAa(Scaffold::AA_O, UniqueName(g, "O")); + sc.SetAa(Scaffold::AA_CB, UniqueName(g, "CB")); + + if (sc.Aa(Scaffold::AA_N) < 0 || sc.Aa(Scaffold::AA_CA) < 0 || + sc.Aa(Scaffold::AA_C) < 0) { + for (int ni = 0; ni < g.Natom(); ni++) { + if (g.Elt(ni) != Atom::NITROGEN) continue; + Iarray const& nNbr = g.HeavyNbr(ni); + for (Iarray::const_iterator ca = nNbr.begin(); ca != nNbr.end(); ++ca) { + if (g.Elt(*ca) != Atom::CARBON) continue; + Iarray const& caNbr = g.HeavyNbr(*ca); + for (Iarray::const_iterator cj = caNbr.begin(); cj != caNbr.end(); ++cj) { + if (*cj == ni || g.Elt(*cj) != Atom::CARBON) continue; + int oIdx = -1; + int nTermO = 0; + Iarray const& cNbr = g.HeavyNbr(*cj); + for (Iarray::const_iterator oj = cNbr.begin(); oj != cNbr.end(); ++oj) { + if (IsTermO(g, *oj)) { nTermO++; oIdx = *oj; } + } + if (nTermO != 1 || (int)cNbr.size() > 3) continue; + if (sc.Aa(Scaffold::AA_N) < 0) sc.SetAa(Scaffold::AA_N, ni); + if (sc.Aa(Scaffold::AA_CA) < 0) sc.SetAa(Scaffold::AA_CA, *ca); + if (sc.Aa(Scaffold::AA_C) < 0) sc.SetAa(Scaffold::AA_C, *cj); + if (sc.Aa(Scaffold::AA_O) < 0) sc.SetAa(Scaffold::AA_O, oIdx); + } + } + } + } + + int n = sc.Aa(Scaffold::AA_N); + int ca = sc.Aa(Scaffold::AA_CA); + int c = sc.Aa(Scaffold::AA_C); + if (sc.Aa(Scaffold::AA_O) < 0 && c >= 0) { + Iarray const& cn = g.HeavyNbr(c); + for (Iarray::const_iterator j = cn.begin(); j != cn.end(); ++j) { + if (IsTermO(g, *j)) { sc.SetAa(Scaffold::AA_O, *j); break; } + } + } + if (sc.Aa(Scaffold::AA_CB) < 0 && ca >= 0) { + Iarray const& can = g.HeavyNbr(ca); + for (Iarray::const_iterator j = can.begin(); j != can.end(); ++j) { + if (g.Elt(*j) == Atom::CARBON && *j != c) + sc.SetAa(Scaffold::AA_CB, *j); + } + } + + n = sc.Aa(Scaffold::AA_N); + ca = sc.Aa(Scaffold::AA_CA); + c = sc.Aa(Scaffold::AA_C); + sc.aaOk_ = (n >= 0 && ca >= 0 && c >= 0); + if (!sc.aaOk_) return; + + sc.isPro_ = false; + if (n >= 0) { + int nC = 0; + Iarray const& nn = g.HeavyNbr(n); + for (Iarray::const_iterator j = nn.begin(); j != nn.end(); ++j) { + if (g.Elt(*j) == Atom::CARBON) nC++; + } + sc.isPro_ = (nC >= 2); + } + + if (sc.kind_ == "unknown") { + sc.kind_ = "amino"; + if (sc.isPro_) sc.family_ = "pro"; + else if (sc.Aa(Scaffold::AA_CB) < 0) sc.family_ = "gly"; + else sc.family_ = "std"; + } +} + +/// Detect nucleotide / sugar / base / amino / unknown and fill Scaffold roles. +/** Order of decisions: + * 1. Find 5- and 6-rings; pick a furanose if present → AssignSugarRoles. + * 2. Glycosidic N on C1′ → base component + purine/pyrimidine family. + * 3. kind_ = nucleotide | sugar | base | unknown from what was found. + * 4. NameHints / ModXNA metadata fill any still-empty roles. + * 5. Guard: clear a false C5′ that is actually exo on C1′ (fragments). + * 6. If still unknown, DetectAmino (peptide motif → kind amino). + * + * Amino acids without a furanose get kind "amino". Ligands (benzene, phenol) + * stay "unknown" and rely on unique-name seeding in MapGraphs. + */ +void DetectScaffold(Graph const& g, Scaffold& sc) { + sc = Scaffold(); + if (g.Natom() == 0) { + sc.notes_.push_back("empty topology"); + return; + } + Cycles rings5 = CyclesOfLength(g, 5); + Cycles rings6 = CyclesOfLength(g, 6); + + Iarray pAtoms; + for (int i = 0; i < g.Natom(); i++) { + if (g.Elt(i) == Atom::PHOSPHORUS) pAtoms.push_back(i); + } + if (!pAtoms.empty()) { + sc.Set(Scaffold::P, pAtoms[0]); + sc.hasP_ = true; + if (pAtoms.size() > 1) + sc.notes_.push_back("multiple P atoms; using first"); + } + + Iarray sugar = PickSugarRing(g, rings5); + if (!sugar.empty()) { + sc.sugarRing_ = sugar; + AssignSugarRoles(g, sc, sugar); + } else { + sc.notes_.push_back("no furanose ring"); + } + + int c1 = sc.Get(Scaffold::C1p); + if (c1 >= 0) { + Iarray nchi; + Iarray const& hn = g.HeavyNbr(c1); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (g.Elt(*j) == Atom::NITROGEN) nchi.push_back(*j); + } + if (!nchi.empty()) { + sc.Set(Scaffold::CHI, nchi[0]); + sc.base_ = BaseAtoms(g, nchi[0], c1); + sc.family_ = BaseFamily(g, nchi[0], rings5, rings6); + } + } + + if (sc.Get(Scaffold::C2p) >= 0) + sc.twoPrime_ = TwoPrimePattern(g, sc); + + if (sc.Get(Scaffold::CHI) >= 0 && !sugar.empty()) { + sc.kind_ = "nucleotide"; + sc.ok_ = true; + } else if (!sugar.empty()) { + sc.kind_ = "sugar"; + sc.ok_ = (sc.Get(Scaffold::C1p) >= 0 && sc.Get(Scaffold::O4p) >= 0); + } else { + int nN = 0, nC = 0; + for (int i = 0; i < g.Natom(); i++) { + if (g.Elt(i) == Atom::NITROGEN) nN++; + if (g.Elt(i) == Atom::CARBON) nC++; + } + bool looksBase = (nN >= 2 && nC >= 4 && (!rings5.empty() || !rings6.empty())); + if (sc.Get(Scaffold::CHI) >= 0 || looksBase) { + sc.kind_ = "base"; + if (sc.Get(Scaffold::CHI) < 0) { + int chi = GuessGlycosidicN(g, rings5, rings6); + sc.Set(Scaffold::CHI, chi); + if (chi >= 0) { + sc.base_ = BaseAtoms(g, chi, -1); + sc.family_ = BaseFamily(g, chi, rings5, rings6); + } + } + sc.ok_ = (sc.Get(Scaffold::CHI) >= 0); + } else { + sc.kind_ = "unknown"; + } + } + + NameHints(g, sc); + ApplyModxna(g, sc); + + int chi = sc.Get(Scaffold::CHI); + if (chi >= 0) { + sc.base_ = BaseAtoms(g, chi, sc.Get(Scaffold::C1p)); + sc.family_ = BaseFamily(g, chi, rings5, rings6); + } + if (sc.Get(Scaffold::C2p) >= 0 && sc.kind_ != "base") + sc.twoPrime_ = TwoPrimePattern(g, sc); + + // A C1′ methyl/base cap must not be labeled C5′ (common on nucleoside fragments). + int c5 = sc.Get(Scaffold::C5p); + c1 = sc.Get(Scaffold::C1p); + int c4 = sc.Get(Scaffold::C4p); + if (c5 >= 0 && c1 >= 0) { + Iarray const& c1n = g.HeavyNbr(c1); + bool onC1 = std::find(c1n.begin(), c1n.end(), c5) != c1n.end(); + bool onC4 = false; + if (c4 >= 0) { + Iarray const& c4n = g.HeavyNbr(c4); + onC4 = std::find(c4n.begin(), c4n.end(), c5) != c4n.end(); + } + if (onC1 && !onC4) { + sc.Set(Scaffold::C5p, -1); + sc.notes_.push_back("cleared C5' attached to C1'"); + c5 = -1; + } + } + + if (sc.kind_ == "unknown") + DetectAmino(g, sc); +} + +// ============================================================================= +// Correspondence — seed, grow, names, hydrogens +// ============================================================================= + +/// Neighbor signature used to grow the map: (element, heavy degree, sorted neighbor elements). +/** S vs Se (Cys/Sec) and O vs H (phenol/benzene at C1) fail this test, which + * is how those substitutions become insertions rather than forced matches. + * Hydrogens are excluded from degree and from the neighbor-element list so + * protonation differences do not block heavy-atom growth. + */ +struct AtomSig { + Atom::AtomicElementType elt_; + int deg_; + std::vector nbr_; + bool operator<(AtomSig const& rhs) const { + if (elt_ != rhs.elt_) return elt_ < rhs.elt_; + if (deg_ != rhs.deg_) return deg_ < rhs.deg_; + return nbr_ < rhs.nbr_; + } +}; + +/// Build the grow signature of atom i (hydrogens are excluded from degree / neighbors). +static AtomSig MakeSig(Graph const& g, int i) { + AtomSig s; + s.elt_ = g.Elt(i); + s.deg_ = g.HeavyDeg(i); + Iarray const& hn = g.HeavyNbr(i); + s.nbr_.reserve(hn.size()); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) + s.nbr_.push_back(g.Elt(*j)); + std::sort(s.nbr_.begin(), s.nbr_.end()); + return s; +} + +/// Stable order for same-signature twins: unique names if possible, else atom index. +static Iarray OrderTwins(Graph const& g, Iarray kids) { + std::set names; + for (Iarray::const_iterator j = kids.begin(); j != kids.end(); ++j) + names.insert(g.Name(*j)); + if ((int)names.size() == (int)kids.size()) { + std::sort(kids.begin(), kids.end(), [&g](int a, int b){ return g.Name(a) < g.Name(b); }); + return kids; + } + std::sort(kids.begin(), kids.end()); + return kids; +} + +/// Record tgt[n] ↔ tpl[r] unless either side is already claimed. Returns true if newly added. +static bool AddMap(Iarray& mapping, Iarray& usedTpl, int n, int r) { + if (n < 0 || r < 0) return false; + if (mapping[n] >= 0) return (mapping[n] == r); + if (usedTpl[r] != 0) return false; + mapping[n] = r; + usedTpl[r] = 1; + return true; +} + +/// Grow the correspondence from already-mapped atoms by unique neighbor signatures. +/** Repeated until a pass adds nothing. For each already-mapped pair (n↔r): + * - Collect unmapped neighbors of n and unused neighbors of r. + * - Bucket by AtomSig. + * - If a signature has exactly one atom on each side → map them. + * - If equal-sized twin groups (H5′/H5″, OP1/OP2) → OrderTwins by name. + * Signatures that appear on only one side are left unmatched (insertions / + * unmatched template). This is how a 2′-OH sits next to C2′ without claiming + * template H2″. + */ +void Grow(Graph const& tgt, Graph const& tpl, + Iarray& mapping, Iarray& usedTpl) +{ + bool changed = true; + while (changed) { + changed = false; + Iarray mappedTgt; + for (int n = 0; n < tgt.Natom(); n++) { + if (mapping[n] >= 0) mappedTgt.push_back(n); + } + for (Iarray::const_iterator it = mappedTgt.begin(); it != mappedTgt.end(); ++it) { + int n = *it; + int r = mapping[n]; + Iarray nUn, rUn; + Iarray const& nn = tgt.Nbr(n); + for (Iarray::const_iterator j = nn.begin(); j != nn.end(); ++j) { + if (mapping[*j] < 0) nUn.push_back(*j); + } + Iarray const& rn = tpl.Nbr(r); + for (Iarray::const_iterator j = rn.begin(); j != rn.end(); ++j) { + if (!usedTpl[*j]) rUn.push_back(*j); + } + if (nUn.empty() || rUn.empty()) continue; + std::map nBy, rBy; + for (Iarray::const_iterator j = nUn.begin(); j != nUn.end(); ++j) + nBy[MakeSig(tgt, *j)].push_back(*j); + for (Iarray::const_iterator j = rUn.begin(); j != rUn.end(); ++j) + rBy[MakeSig(tpl, *j)].push_back(*j); + for (std::map::iterator nb = nBy.begin(); nb != nBy.end(); ++nb) { + std::map::iterator rb = rBy.find(nb->first); + if (rb == rBy.end()) continue; + Iarray& nlist = nb->second; + Iarray& rlist = rb->second; + if (nlist.size() == 1 && rlist.size() == 1) { + if (AddMap(mapping, usedTpl, nlist[0], rlist[0])) + changed = true; + } else if (nlist.size() == rlist.size() && nlist.size() > 1) { + Iarray ns = OrderTwins(tgt, nlist); + Iarray rs = OrderTwins(tpl, rlist); + for (size_t k = 0; k < ns.size(); k++) { + if (AddMap(mapping, usedTpl, ns[k], rs[k])) + changed = true; + } + } + } + } + } +} + +/// Pair leftover atoms that share a unique name *and* the same element. +/** Element check is required: mol2 atom name SE would otherwise look like + * sulfur, and even with OFF, SE must not be forced onto SG. + */ +void MatchUniqueNames(Graph const& tgt, Graph const& tpl, + Iarray& mapping, Iarray& usedTpl) +{ + std::map refByName; + for (int i = 0; i < tpl.Natom(); i++) + refByName[tpl.Name(i)].push_back(i); + for (int i = 0; i < tgt.Natom(); i++) { + if (mapping[i] >= 0) continue; + std::map::const_iterator it = refByName.find(tgt.Name(i)); + if (it == refByName.end()) continue; + Iarray hits; + for (Iarray::const_iterator r = it->second.begin(); r != it->second.end(); ++r) { + if (!usedTpl[*r]) hits.push_back(*r); + } + if (hits.size() == 1 && tgt.Elt(i) == tpl.Elt(hits[0])) + AddMap(mapping, usedTpl, i, hits[0]); + } +} + +/// Zip hydrogens of already-mapped heavies, ordered by name (H61 before H62). +void MatchHydrogens(Graph const& tgt, Graph const& tpl, + Iarray& mapping, Iarray& usedTpl) +{ + for (int n = 0; n < tgt.Natom(); n++) { + int r = mapping[n]; + if (r < 0) continue; + Iarray nH, rH; + Iarray const& th = tgt.HNbr(n); + for (Iarray::const_iterator j = th.begin(); j != th.end(); ++j) { + if (mapping[*j] < 0) nH.push_back(*j); + } + Iarray const& rh = tpl.HNbr(r); + for (Iarray::const_iterator j = rh.begin(); j != rh.end(); ++j) { + if (!usedTpl[*j]) rH.push_back(*j); + } + if (nH.empty() || rH.empty()) continue; + Iarray ns = OrderTwins(tgt, nH); + Iarray rs = OrderTwins(tpl, rH); + size_t nzip = std::min(ns.size(), rs.size()); + for (size_t k = 0; k < nzip; k++) + AddMap(mapping, usedTpl, ns[k], rs[k]); + } +} + +/// Seed, grow, unique-name, grow, hydrogens. mapping[tgt] = tpl or −1. +/** SEED_AUTO picking order: + * 1. SEED_NA if either scaffold looks nucleotide-like (ok_, P, C1′, or χ). + * 2. Else SEED_AA if either has aaOk_ (peptide backbone). + * 3. Else SEED_NAMES (unique atom names only). + * SEED_NONE skips all pairing (debug / empty map). + * + * After the seed pass we always: Grow → MatchUniqueNames → Grow → MatchHydrogens + * (unless SEED_NONE). The second Grow catches atoms that only become unique + * after names filled gaps; hydrogens run last so heavy topology is fixed. + * + * \param mapping resized to tgt.Natom(); −1 = insertion. + */ +void MapGraphs(Graph const& tgt, Graph const& tpl, + Scaffold const& tgtSc, Scaffold const& tplSc, + Iarray& mapping, TIMatch::SeedType seedIn) +{ + mapping.assign(tgt.Natom(), -1); + Iarray usedTpl(tpl.Natom(), 0); + + TIMatch::SeedType mode = seedIn; + if (mode == TIMatch::SEED_AUTO) { + bool na = tgtSc.ok_ || tplSc.ok_ || + tgtSc.Get(Scaffold::P) >= 0 || tplSc.Get(Scaffold::P) >= 0 || + tgtSc.Get(Scaffold::C1p) >= 0 || tplSc.Get(Scaffold::C1p) >= 0 || + tgtSc.Get(Scaffold::CHI) >= 0 || tplSc.Get(Scaffold::CHI) >= 0; + bool aa = tgtSc.aaOk_ || tplSc.aaOk_; + if (na) mode = TIMatch::SEED_NA; + else if (aa) mode = TIMatch::SEED_AA; + else mode = TIMatch::SEED_NAMES; + } + + if (mode == TIMatch::SEED_NA) { + static const Scaffold::Role SEED[] = { + Scaffold::P, Scaffold::O5p, Scaffold::C5p, Scaffold::C4p, Scaffold::O4p, + Scaffold::C1p, Scaffold::C3p, Scaffold::C2p, Scaffold::O3p, Scaffold::CHI + }; + for (int k = 0; k < 10; k++) + AddMap(mapping, usedTpl, tgtSc.Get(SEED[k]), tplSc.Get(SEED[k])); + size_t nop = std::min(tgtSc.op_.size(), tplSc.op_.size()); + for (size_t i = 0; i < nop; i++) + AddMap(mapping, usedTpl, tgtSc.op_[i], tplSc.op_[i]); + } + + if (mode == TIMatch::SEED_AA) { + static const Scaffold::AaRole SEED[] = { + Scaffold::AA_N, Scaffold::AA_CA, Scaffold::AA_C, Scaffold::AA_O, Scaffold::AA_CB + }; + for (int k = 0; k < 5; k++) + AddMap(mapping, usedTpl, tgtSc.Aa(SEED[k]), tplSc.Aa(SEED[k])); + } + + if (mode != TIMatch::SEED_NONE) { + Grow(tgt, tpl, mapping, usedTpl); + MatchUniqueNames(tgt, tpl, mapping, usedTpl); + Grow(tgt, tpl, mapping, usedTpl); + MatchHydrogens(tgt, tpl, mapping, usedTpl); + } +} + +// ============================================================================= +// Output layouts — target permutation and dual-topology slots +// ============================================================================= + +/// Depth-first emit of an unmapped insertion subtree (heavy atoms before hydrogens). +/** used[] marks atoms already placed in \p order. mappedNew[] marks target + * atoms that correspond to a template atom — those are *not* walked here + * (they appear when OutputOrder / DualOrder visits their template parent). + */ +static void EmitTree(Graph const& g, int n, int parent, + std::vector& used, Iarray const& mappedNew, + Iarray& order) +{ + if (used[n]) return; + used[n] = 1; + order.push_back(n); + Iarray kids; + Iarray const& nbr = g.Nbr(n); + for (Iarray::const_iterator j = nbr.begin(); j != nbr.end(); ++j) { + if (*j == parent || used[*j] || mappedNew[*j]) continue; + kids.push_back(*j); + } + std::sort(kids.begin(), kids.end(), [&g](int a, int b) { + if (g.IsHydrogen(a) != g.IsHydrogen(b)) return !g.IsHydrogen(a) && g.IsHydrogen(b); + if (g.Elt(a) != g.Elt(b)) return g.Elt(a) < g.Elt(b); + return g.Name(a) < g.Name(b); + }); + for (Iarray::const_iterator j = kids.begin(); j != kids.end(); ++j) + EmitTree(g, *j, n, used, mappedNew, order); +} + +/// Permutation of target atoms: walk the template, emit the mapped atom, then insertions. +/** Algorithm: + * For each template atom r in parentOrder (file / NA / AA walk): + * If some target n maps to r, append n, then EmitTree any unmapped + * neighbors of n (insertions bonded to that parent — O2′ after C2′, + * phenolic OH after C1, Se after CB). + * Unmatched template atoms are *skipped* (H2″ of dA, H1 of benzene, SG). + * Remaining unused target atoms (unattached leftovers) insert just before + * the O3′ / anchor slot, or at the end if no anchor. + * + * Return value is Map[newIndex] = old target atom (ModifyByMap convention). + */ +Iarray OutputOrder(Graph const& tgt, Graph const& tpl, + Scaffold const& tplSc, + Iarray const& mapping, + Iarray const& parentOrder, + std::string const& anchorName) +{ + Iarray refToNew(tpl.Natom(), -1); + Iarray mappedNew(tgt.Natom(), 0); + for (int n = 0; n < tgt.Natom(); n++) { + if (mapping[n] >= 0) { + refToNew[mapping[n]] = n; + mappedNew[n] = 1; + } + } + std::vector used(tgt.Natom(), 0); + Iarray order; + int o3slot = -1; + int o3 = tplSc.Get(Scaffold::O3p); + if (o3 < 0 && !anchorName.empty()) { + Iarray hits = tpl.FindByName(anchorName); + if (hits.size() == 1) o3 = hits[0]; + } + + for (Iarray::const_iterator pr = parentOrder.begin(); pr != parentOrder.end(); ++pr) { + int r = *pr; + if (r == o3) o3slot = (int)order.size(); + if (r < 0 || r >= tpl.Natom() || refToNew[r] < 0) continue; + int n = refToNew[r]; + if (used[n]) continue; + used[n] = 1; + order.push_back(n); + Iarray kids; + Iarray const& nbr = tgt.Nbr(n); + for (Iarray::const_iterator j = nbr.begin(); j != nbr.end(); ++j) { + if (!used[*j] && !mappedNew[*j]) kids.push_back(*j); + } + std::sort(kids.begin(), kids.end(), [&tgt](int a, int b) { + if (tgt.IsHydrogen(a) != tgt.IsHydrogen(b)) + return !tgt.IsHydrogen(a) && tgt.IsHydrogen(b); + if (tgt.Elt(a) != tgt.Elt(b)) return tgt.Elt(a) < tgt.Elt(b); + return tgt.Name(a) < tgt.Name(b); + }); + for (Iarray::const_iterator j = kids.begin(); j != kids.end(); ++j) + EmitTree(tgt, *j, n, used, mappedNew, order); + } + + Iarray leftover; + for (int i = 0; i < tgt.Natom(); i++) { + if (!used[i]) leftover.push_back(i); + } + std::sort(leftover.begin(), leftover.end(), [&tgt](int a, int b) { + if (tgt.IsHydrogen(a) != tgt.IsHydrogen(b)) + return !tgt.IsHydrogen(a) && tgt.IsHydrogen(b); + if (tgt.Elt(a) != tgt.Elt(b)) return tgt.Elt(a) < tgt.Elt(b); + if (tgt.Name(a) != tgt.Name(b)) return tgt.Name(a) < tgt.Name(b); + return a < b; + }); + if (!leftover.empty()) { + int insertAt = (o3slot >= 0) ? o3slot : (int)order.size(); + order.insert(order.begin() + insertAt, leftover.begin(), leftover.end()); + } + return order; +} + +/// Dual-topology layout: every template atom keeps a slot; insertions are extra slots. +/** Same walk as OutputOrder, but when template atom r has no target partner we + * still emit DualSlot(r, −1) (real on λ=0 / dummy on λ=1). Target-only + * insertions are DualSlot(−1, n) (dummy on λ=0 / real on λ=1). Shared atoms + * are DualSlot(r, n). + * + * Size invariant: dual.size() == nMapped + nInsertion + nUnmappedTpl. + * Exec_TIMap::BuildTiUnit consumes this vector to build equal-NATOM end states. + */ +static std::vector +DualOrder(Graph const& tgt, Graph const& tpl, + Scaffold const& tplSc, + Iarray const& mapping, + Iarray const& parentOrder, + std::string const& anchorName) +{ + typedef TIMatch::Result::DualSlot Slot; + Iarray refToNew(tpl.Natom(), -1); + Iarray mappedNew(tgt.Natom(), 0); + for (int n = 0; n < tgt.Natom(); n++) { + if (mapping[n] >= 0) { + refToNew[mapping[n]] = n; + mappedNew[n] = 1; + } + } + std::vector used(tgt.Natom(), 0); + std::vector dual; + int o3slot = -1; + int o3 = tplSc.Get(Scaffold::O3p); + if (o3 < 0 && !anchorName.empty()) { + Iarray hits = tpl.FindByName(anchorName); + if (hits.size() == 1) o3 = hits[0]; + } + + for (Iarray::const_iterator pr = parentOrder.begin(); pr != parentOrder.end(); ++pr) { + int r = *pr; + if (r == o3) o3slot = (int)dual.size(); + if (r < 0 || r >= tpl.Natom()) continue; + if (refToNew[r] < 0) { + dual.push_back(Slot(r, -1)); + continue; + } + int n = refToNew[r]; + if (used[n]) continue; + used[n] = 1; + dual.push_back(Slot(r, n)); + Iarray kids; + Iarray const& nbr = tgt.Nbr(n); + for (Iarray::const_iterator j = nbr.begin(); j != nbr.end(); ++j) { + if (!used[*j] && !mappedNew[*j]) kids.push_back(*j); + } + std::sort(kids.begin(), kids.end(), [&tgt](int a, int b) { + if (tgt.IsHydrogen(a) != tgt.IsHydrogen(b)) + return !tgt.IsHydrogen(a) && tgt.IsHydrogen(b); + if (tgt.Elt(a) != tgt.Elt(b)) return tgt.Elt(a) < tgt.Elt(b); + return tgt.Name(a) < tgt.Name(b); + }); + Iarray tree; + for (Iarray::const_iterator j = kids.begin(); j != kids.end(); ++j) + EmitTree(tgt, *j, n, used, mappedNew, tree); + for (Iarray::const_iterator j = tree.begin(); j != tree.end(); ++j) + dual.push_back(Slot(-1, *j)); + } + + Iarray leftover; + for (int i = 0; i < tgt.Natom(); i++) { + if (!used[i]) leftover.push_back(i); + } + std::sort(leftover.begin(), leftover.end(), [&tgt](int a, int b) { + if (tgt.IsHydrogen(a) != tgt.IsHydrogen(b)) + return !tgt.IsHydrogen(a) && tgt.IsHydrogen(b); + if (tgt.Elt(a) != tgt.Elt(b)) return tgt.Elt(a) < tgt.Elt(b); + if (tgt.Name(a) != tgt.Name(b)) return tgt.Name(a) < tgt.Name(b); + return a < b; + }); + if (!leftover.empty()) { + int insertAt = (o3slot >= 0) ? o3slot : (int)dual.size(); + std::vector extra; + extra.reserve(leftover.size()); + for (Iarray::const_iterator j = leftover.begin(); j != leftover.end(); ++j) + extra.push_back(Slot(-1, *j)); + dual.insert(dual.begin() + insertAt, extra.begin(), extra.end()); + } + return dual; +} + +// ============================================================================= +// Canonical walks — invent Amber-like order when no parent file order is used +// ============================================================================= + +/// Mark atom i used and append it; no-op if i < 0 or already emitted. +static void EmitIdx(int i, std::vector& used, Iarray& order) { + if (i < 0 || used[i]) return; + used[i] = 1; + order.push_back(i); +} + +/// Emit a heavy atom and its hydrogens (name-sorted) for the canonical NA walk. +static void EmitHeavyAndH(Graph const& g, int i, + std::vector& used, Iarray& order) +{ + EmitIdx(i, used, order); + if (i < 0) return; + Iarray hs = g.HNbr(i); + std::sort(hs.begin(), hs.end(), [&g](int a, int b){ return g.Name(a) < g.Name(b); }); + for (Iarray::const_iterator h = hs.begin(); h != hs.end(); ++h) + EmitIdx(*h, used, order); +} + +/// Lower number = earlier in the nucleobase BFS (ring atoms before exo O/N/C, H last). +static int BasePriority(Graph const& g, int i) { + if (g.IsHydrogen(i)) return 9; + int nRingish = g.HeavyDeg(i); + if (g.Elt(i) == Atom::OXYGEN && nRingish <= 1) return 5; + if (g.Elt(i) == Atom::NITROGEN && nRingish <= 1) return 4; + if (g.Elt(i) == Atom::CARBON && nRingish <= 1) return 6; + return 0; +} + +/// BFS of the nucleobase from χ, never crossing back onto C1′. +static Iarray WalkBase(Graph const& g, int chi, int c1) { + std::set skip; + if (c1 >= 0) skip.insert(c1); + std::set seen; + seen.insert(chi); + Iarray order(1, chi); + Iarray stack(1, chi); + size_t q = 0; + while (q < stack.size()) { + int i = stack[q++]; + Iarray kids; + Iarray const& hn = g.HeavyNbr(i); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (seen.find(*j) != seen.end()) continue; + if (skip.find(*j) != skip.end()) continue; + kids.push_back(*j); + } + std::sort(kids.begin(), kids.end(), [&g](int a, int b) { + int pa = BasePriority(g, a); + int pb = BasePriority(g, b); + if (pa != pb) return pa < pb; + if (g.Name(a) != g.Name(b)) return g.Name(a) < g.Name(b); + return a < b; + }); + for (Iarray::const_iterator j = kids.begin(); j != kids.end(); ++j) { + seen.insert(*j); + order.push_back(*j); + stack.push_back(*j); + } + } + return order; +} + +/// Emit atom i and every still-unused neighbor subtree (used for 2′ substituents). +static void EmitSubtree(Graph const& g, int i, int parent, + std::vector& used, Iarray& order) +{ + if (used[i]) return; + used[i] = 1; + order.push_back(i); + Iarray kids; + Iarray const& nbr = g.Nbr(i); + for (Iarray::const_iterator j = nbr.begin(); j != nbr.end(); ++j) { + if (*j != parent && !used[*j]) kids.push_back(*j); + } + std::sort(kids.begin(), kids.end(), [&g](int a, int b) { + if (g.IsHydrogen(a) != g.IsHydrogen(b)) + return !g.IsHydrogen(a) && g.IsHydrogen(b); + if (g.Elt(a) != g.Elt(b)) return g.Elt(a) < g.Elt(b); + return g.Name(a) < g.Name(b); + }); + for (Iarray::const_iterator j = kids.begin(); j != kids.end(); ++j) + EmitSubtree(g, *j, i, used, order); +} + +/// Canonical nucleic-acid walk used when the user requests naorder. +/** Intended shared-atom order for ModXNA-class / Amber nucleotides: + * + * P → OP (+H, name-sorted) → O5′ → C5′ → C4′ → O4′ → C1′ → + * base BFS from χ (never recrossing C1′) → + * C3′ → C2′ → 2′ substituent subtrees → O3′ → leftovers + * + * Hydrogens follow their heavy atom (name-sorted). Leftovers append in + * index order. See manual section "How the nucleic-acid walk was chosen". + */ +Iarray CanonicalWalk(Graph const& g, Scaffold const& sc) { + std::vector used(g.Natom(), 0); + Iarray order; + EmitIdx(sc.Get(Scaffold::P), used, order); + for (Iarray::const_iterator op = sc.op_.begin(); op != sc.op_.end(); ++op) + EmitHeavyAndH(g, *op, used, order); + EmitHeavyAndH(g, sc.Get(Scaffold::O5p), used, order); + EmitHeavyAndH(g, sc.Get(Scaffold::C5p), used, order); + EmitHeavyAndH(g, sc.Get(Scaffold::C4p), used, order); + EmitHeavyAndH(g, sc.Get(Scaffold::O4p), used, order); + EmitHeavyAndH(g, sc.Get(Scaffold::C1p), used, order); + int chi = sc.Get(Scaffold::CHI); + if (chi >= 0) { + Iarray base = WalkBase(g, chi, sc.Get(Scaffold::C1p)); + for (Iarray::const_iterator i = base.begin(); i != base.end(); ++i) + EmitHeavyAndH(g, *i, used, order); + } + EmitHeavyAndH(g, sc.Get(Scaffold::C3p), used, order); + int c2 = sc.Get(Scaffold::C2p); + EmitHeavyAndH(g, c2, used, order); + if (c2 >= 0) { + std::set skip; + if (sc.Get(Scaffold::C1p) >= 0) skip.insert(sc.Get(Scaffold::C1p)); + if (sc.Get(Scaffold::C3p) >= 0) skip.insert(sc.Get(Scaffold::C3p)); + if (sc.Get(Scaffold::O4p) >= 0) skip.insert(sc.Get(Scaffold::O4p)); + Iarray kids; + Iarray const& nbr = g.Nbr(c2); + for (Iarray::const_iterator j = nbr.begin(); j != nbr.end(); ++j) { + if (!used[*j] && skip.find(*j) == skip.end()) kids.push_back(*j); + } + std::sort(kids.begin(), kids.end(), [&g](int a, int b) { + if (g.IsHydrogen(a) != g.IsHydrogen(b)) + return !g.IsHydrogen(a) && g.IsHydrogen(b); + if (g.Elt(a) != g.Elt(b)) return g.Elt(a) < g.Elt(b); + return g.Name(a) < g.Name(b); + }); + for (Iarray::const_iterator j = kids.begin(); j != kids.end(); ++j) + EmitSubtree(g, *j, c2, used, order); + } + EmitHeavyAndH(g, sc.Get(Scaffold::O3p), used, order); + for (int i = 0; i < g.Natom(); i++) { + if (!used[i]) EmitIdx(i, used, order); + } + return order; +} + +/// Count unused heavy atoms reachable from start without crossing parent or used. +static int HeavySubtreeSize(Graph const& g, int start, int parent, + std::vector const& used) +{ + std::vector seen(g.Natom(), 0); + Iarray stack(1, start); + seen[start] = 1; + int count = 0; + while (!stack.empty()) { + int i = stack.back(); + stack.pop_back(); + if (used[i] && i != start) continue; + if (!g.IsHydrogen(i)) count++; + Iarray const& nbr = g.Nbr(i); + for (Iarray::const_iterator j = nbr.begin(); j != nbr.end(); ++j) { + if (*j == parent || seen[*j] || (used[*j] && *j != start)) continue; + seen[*j] = 1; + stack.push_back(*j); + } + } + return count; +} + +/// Walk a side-chain branch the way amino19.lib does: heavy, its hydrogens, then +/// child heavies (smaller unused subtree first, then name). +static void EmitAaBranch(Graph const& g, int i, int parent, + std::vector& used, Iarray& order, + std::set const& skip) +{ + if (i < 0 || used[i] || skip.find(i) != skip.end()) return; + EmitHeavyAndH(g, i, used, order); + Iarray kids; + Iarray const& hn = g.HeavyNbr(i); + for (Iarray::const_iterator j = hn.begin(); j != hn.end(); ++j) { + if (*j == parent || used[*j] || skip.find(*j) != skip.end()) continue; + kids.push_back(*j); + } + std::sort(kids.begin(), kids.end(), [&g, &used, i](int a, int b) { + int sa = HeavySubtreeSize(g, a, i, used); + int sb = HeavySubtreeSize(g, b, i, used); + if (sa != sb) return sa < sb; + return g.Name(a) < g.Name(b); + }); + for (Iarray::const_iterator j = kids.begin(); j != kids.end(); ++j) + EmitAaBranch(g, *j, i, used, order, skip); +} + +/// Canonical amino-acid walk used when the user requests aaorder (ff19SB). +/** N (+ amide H) → [Pro: CD-ring to CB] → CA (+ HA) → CB side chain → C → O. + * Leftovers (OXT, caps) append in index order. + */ +Iarray CanonicalAaWalk(Graph const& g, Scaffold const& sc) { + std::vector used(g.Natom(), 0); + Iarray order; + int n = sc.Aa(Scaffold::AA_N); + int ca = sc.Aa(Scaffold::AA_CA); + int c = sc.Aa(Scaffold::AA_C); + int o = sc.Aa(Scaffold::AA_O); + int cb = sc.Aa(Scaffold::AA_CB); + std::set skip; + if (n >= 0) skip.insert(n); + if (ca >= 0) skip.insert(ca); + if (c >= 0) skip.insert(c); + if (o >= 0) skip.insert(o); + + EmitHeavyAndH(g, n, used, order); + if (sc.isPro_ && n >= 0) { + Iarray const& nn = g.HeavyNbr(n); + for (Iarray::const_iterator j = nn.begin(); j != nn.end(); ++j) { + if (*j != ca && !used[*j]) + EmitAaBranch(g, *j, n, used, order, skip); + } + } + EmitHeavyAndH(g, ca, used, order); + if (!sc.isPro_) { + if (cb >= 0) + EmitAaBranch(g, cb, ca, used, order, skip); + else if (ca >= 0) { + Iarray const& can = g.HeavyNbr(ca); + for (Iarray::const_iterator j = can.begin(); j != can.end(); ++j) { + if (*j != n && *j != c && !used[*j] && !g.IsHydrogen(*j)) + EmitAaBranch(g, *j, ca, used, order, skip); + } + } + } + EmitHeavyAndH(g, c, used, order); + EmitHeavyAndH(g, o, used, order); + for (int i = 0; i < g.Natom(); i++) { + if (!used[i]) EmitIdx(i, used, order); + } + return order; +} + +/// Template order: file order, NA walk, or amino-acid walk (ff19SB). +/** Called once per Match() after MapGraphs. The returned indices are template + * atom indices in the order OutputOrder / DualOrder should visit them. + * useNaOrder / useAaOrder come from timap naorder / aaorder (mutually exclusive + * at the Exec_TIMap layer). If both false, identity 0..N−1 = file order. + */ +Iarray ParentOrder(Graph const& tpl, Scaffold const& sc, bool useNaOrder, bool useAaOrder) { + if (useNaOrder) + return CanonicalWalk(tpl, sc); + if (useAaOrder) + return CanonicalAaWalk(tpl, sc); + Iarray order(tpl.Natom()); + for (int i = 0; i < tpl.Natom(); i++) + order[i] = i; + return order; +} + +} // namespace + +// ============================================================================= +// TIMatch public API +// ============================================================================= + +/// Public wrapper: detect the scaffold of \p top and return CanonicalWalk. +int TIMatch::CanonicalNaOrder(Topology const& top, Iarray& order) const { + Graph g(top); + Scaffold sc; + DetectScaffold(g, sc); + order = CanonicalWalk(g, sc); + return 0; +} + +/// Public wrapper: detect peptide scaffold and return CanonicalAaWalk. +int TIMatch::CanonicalAaOrder(Topology const& top, Iarray& order) const { + Graph g(top); + Scaffold sc; + DetectScaffold(g, sc); + order = CanonicalAaWalk(g, sc); + return 0; +} + +/// Detect scaffolds, map graphs, build the output permutation, count partial-map stats. +/** Steps: + * 1. Reject empty topologies. + * 2. DetectScaffold on tgt and tpl (kinds + notes go into Result). + * 3. MapGraphs → mapping_ (partial correspondence). + * 4. ParentOrder → OutputOrder → DualOrder. + * 5. Assert outputOrder_ is a permutation; dual_ size matches counts. + * 6. Fill nMapped_ / nInsertion_ / nUnmappedTpl_. + * + * Return 0 on success (partial maps are success). Nonzero only on empty + * input or internal invariant failure. + */ +int TIMatch::Match(Topology const& tgtTop, Topology const& tplTop, Result& out) const { + out = Result(); + if (tgtTop.Natom() < 1 || tplTop.Natom() < 1) { + mprinterr("Error: timap: empty topology.\n"); + return 1; + } + Graph tgt(tgtTop); + Graph tpl(tplTop); + Scaffold tgtSc, tplSc; + DetectScaffold(tgt, tgtSc); + DetectScaffold(tpl, tplSc); + out.tgtKind_ = tgtSc.kind_; + out.tplKind_ = tplSc.kind_; + out.notes_.insert(out.notes_.end(), tgtSc.notes_.begin(), tgtSc.notes_.end()); + + if (debug_ > 0) { + mprintf("DEBUG: timap tgt kind=%s family=%s 2'=%s ok=%i\n", + tgtSc.kind_.c_str(), tgtSc.family_.c_str(), tgtSc.twoPrime_.c_str(), (int)tgtSc.ok_); + mprintf("DEBUG: timap tpl kind=%s family=%s 2'=%s ok=%i\n", + tplSc.kind_.c_str(), tplSc.family_.c_str(), tplSc.twoPrime_.c_str(), (int)tplSc.ok_); + for (int r = 0; r < Scaffold::NROLES; r++) { + int ti = tgtSc.Get((Scaffold::Role)r); + int pi = tplSc.Get((Scaffold::Role)r); + mprintf("DEBUG: %-4s tgt=%s tpl=%s\n", Scaffold::RoleStr((Scaffold::Role)r), + ti < 0 ? "-" : tgt.Name(ti).c_str(), + pi < 0 ? "-" : tpl.Name(pi).c_str()); + } + } + + MapGraphs(tgt, tpl, tgtSc, tplSc, out.mapping_, seed_); + // Template visit order for OutputOrder / DualOrder (file, NA walk, or AA walk). + Iarray parent = ParentOrder(tpl, tplSc, useNaOrder_, useAaOrder_); + out.outputOrder_ = OutputOrder(tgt, tpl, tplSc, out.mapping_, parent, anchorName_); + out.dual_ = DualOrder(tgt, tpl, tplSc, out.mapping_, parent, anchorName_); + + // ModifyByMap requires a true permutation of the target. + if ((int)out.outputOrder_.size() != tgt.Natom()) { + mprinterr("Error: timap: output order size %zu != %i atoms.\n", + out.outputOrder_.size(), tgt.Natom()); + return 1; + } + std::vector seen(tgt.Natom(), 0); + for (Iarray::const_iterator it = out.outputOrder_.begin(); it != out.outputOrder_.end(); ++it) { + if (*it < 0 || *it >= tgt.Natom() || seen[*it]) { + mprinterr("Error: timap: output order is not a permutation.\n"); + return 1; + } + seen[*it] = 1; + } + + // Partial-map statistics for map headers and series parent scoring. + Iarray usedTpl(tpl.Natom(), 0); + for (int n = 0; n < tgt.Natom(); n++) { + if (out.mapping_[n] >= 0) { + out.nMapped_++; + usedTpl[out.mapping_[n]] = 1; + } else { + out.nInsertion_++; + } + } + for (int r = 0; r < tpl.Natom(); r++) { + if (!usedTpl[r]) out.nUnmappedTpl_++; + } + // Dual layout must account for every shared, insertion, and unmatched tpl atom. + int expectDual = out.nMapped_ + out.nInsertion_ + out.nUnmappedTpl_; + if ((int)out.dual_.size() != expectDual) { + mprinterr("Error: timap: dual-topology size %zu != %i (mapped+ins+unmapped).\n", + out.dual_.size(), expectDual); + return 1; + } + return 0; +} diff --git a/src/TIMatch.h b/src/TIMatch.h new file mode 100644 index 0000000000..6e4f050e68 --- /dev/null +++ b/src/TIMatch.h @@ -0,0 +1,95 @@ +#ifndef INC_TIMATCH_H +#define INC_TIMATCH_H +#include +#include +class Topology; +/// Map a target topology onto a user-supplied template (engine for the timap command). +/** Source encoding is UTF-8. Unicode (O5′, χ, λ, 2′) appears in comments only; + * string literals in the implementation stay ASCII so the file compiles + * without a compiler UTF-8 flag. + * + * The template atom order *is* the thermodynamic-integration (TI) shared-atom + * order. Target atoms that correspond to a template atom occupy that same + * slot. Extra atoms on the target — a 2′-OH on RNA vs DNA, a phenolic OH, + * selenium in place of sulfur — are insertions. Insertions that are bonded + * to a mapped parent are emitted immediately after that parent. Completely + * unmatched leftovers are placed just before an anchor atom (default O3′). + * + * Matching is graph-based and is expected to be partial (unlike atommap, + * which requires a 1:1 correspondence): + * 1. Seed by unique atom names, nucleic-acid scaffold roles, or amino-acid + * N / CA / C / O / CB (ff19SB amino19.lib baseline). + * 2. Grow by pairing unmatched neighbors that share a unique signature + * (element, heavy degree, sorted neighbor elements). + * 3. Attach hydrogens to already-mapped heavy atoms, ordered by name. + * + * The same engine is used for nucleotides, amino acids, and small molecules. + * Nucleic-acid scaffold detection is opportunistic: if no furanose / χ is + * found, the residue is checked for a peptide N–CA–C=O motif (kind "amino"); + * otherwise kind is "unknown" and unique-name seeding is used. + * + * Implementation details, invariants, and "where to change what" live in the + * file-level comment at the top of TIMatch.cpp - start there when maintaining. + */ +class TIMatch { + public: + typedef std::vector Iarray; + /// How the correspondence is started before the grow pass. + enum SeedType { SEED_AUTO = 0, SEED_NAMES, SEED_NA, SEED_AA, SEED_NONE }; + + /// Outcome of Match(): correspondence plus the permutation of the target. + class Result { + public: + Result() : nMapped_(0), nInsertion_(0), nUnmappedTpl_(0) {} + Iarray mapping_; ///< tgt index → template index; −1 = insertion + Iarray outputOrder_; ///< Map[newatom] = old tgt atom (full permutation) + /// Dual-topology slots: tpl-only (dummy on λ=1), tgt-only (dummy on λ=0), or shared. + struct DualSlot { + DualSlot() : tpl_(-1), tgt_(-1) {} + DualSlot(int r, int n) : tpl_(r), tgt_(n) {} + int tpl_; ///< template atom, or −1 if this slot is a target insertion + int tgt_; ///< target atom, or −1 if this slot is an unmatched template atom + bool IsShared() const { return tpl_ >= 0 && tgt_ >= 0; } + bool IsTplOnly() const { return tpl_ >= 0 && tgt_ < 0; } + bool IsTgtOnly() const { return tpl_ < 0 && tgt_ >= 0; } + }; + std::vector dual_; + int nMapped_; ///< target atoms that correspond to a template atom + int nInsertion_; ///< target atoms with no template partner + int nUnmappedTpl_; ///< template atoms with no target partner + std::string tgtKind_; ///< nucleotide | sugar | base | amino | unknown + std::string tplKind_; + std::vector notes_; ///< scaffold warnings (empty sugar, …) + }; + + TIMatch(); + + void SetSeed(SeedType s) { seed_ = s; } + void SetDebug(int d) { debug_ = d; } + void SetUseNaOrder(bool b) { useNaOrder_ = b; } + void SetUseAaOrder(bool b) { useAaOrder_ = b; } + void SetAnchorName(std::string const& n) { anchorName_ = n; } + + /// Map tgt onto tpl. Partial maps are success (return 0) as long as the + /// output order is a permutation of the target. + int Match(Topology const& tgt, Topology const& tpl, Result& out) const; + /// Nucleic-acid canonical walk of top: P → OP → O5′ → C5′ → C4′ → O4′ → + /// C1′ → base from χ → C3′ → C2′ (+ 2′ substituents) → O3′. + int CanonicalNaOrder(Topology const& top, Iarray& order) const; + /// Amino-acid canonical walk of top (ff19SB amino19.lib): N → H → CA → HA → + /// side chain from CB (each heavy, then its hydrogens) → C → O. + /// Proline walks N → CD → … → CB → CA → C → O. + int CanonicalAaOrder(Topology const& top, Iarray& order) const; + + static const char* SeedStr(SeedType); + static SeedType SeedFromString(std::string const&); + + private: + SeedType seed_; + int debug_; + bool useNaOrder_; + bool useAaOrder_; + std::string anchorName_; +}; +#endif + diff --git a/src/cpptrajfiles b/src/cpptrajfiles index bec4fc7f47..ac4f2fc6e2 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -329,6 +329,7 @@ COMMON_SOURCES= \ Exec_Source.cpp \ Exec_SplitCoords.cpp \ Exec_System.cpp \ + Exec_TIMap.cpp \ Exec_Top.cpp \ Exec_Traj.cpp \ Exec_UpdateParameters.cpp \ @@ -437,6 +438,7 @@ COMMON_SOURCES= \ StringRoutines.cpp \ StructureCheck.cpp \ StructureMapper.cpp \ + TIMatch.cpp \ SymmetricRmsdCalc.cpp \ TextBlockBuffer.cpp \ TextFormat.cpp \ diff --git a/src/cpptrajheaders b/src/cpptrajheaders index 3ed37c33a9..6c9e2295c5 100644 --- a/src/cpptrajheaders +++ b/src/cpptrajheaders @@ -346,6 +346,7 @@ CPPTRAJ_HEADERS = \ ./Exec_Source.h \ ./Exec_SplitCoords.h \ ./Exec_System.h \ + ./Exec_TIMap.h \ ./Exec_Top.h \ ./Exec_Traj.h \ ./Exec_UpdateParameters.h \ @@ -487,6 +488,7 @@ CPPTRAJ_HEADERS = \ ./StringRoutines.h \ ./StructureCheck.h \ ./StructureMapper.h \ + ./TIMatch.h \ ./SymbolExporting.h \ ./SymmetricRmsdCalc.h \ ./SymmetricTensor.h \ diff --git a/test/Makefile b/test/Makefile index c4739b09a9..d4ed89f4ed 100644 --- a/test/Makefile +++ b/test/Makefile @@ -77,6 +77,11 @@ test.pucker: test.atommap: @-cd Test_AtomMap && ./RunTest.sh $(OPT) +test.timap: + @-cd Test_TIMap && ./RunTest.sh $(OPT) + +test.timatch: test.timap + test.2drms: @-cd Test_2DRMS && ./RunTest.sh $(OPT) @@ -610,6 +615,7 @@ COMPLETETESTS=test.general \ test.average \ test.pucker \ test.atommap \ + test.timap \ test.datafile \ test.radial \ test.hist \ diff --git a/test/Test_TIMap/5PO.mol2 b/test/Test_TIMap/5PO.mol2 new file mode 100644 index 0000000000..9c5681c6a5 --- /dev/null +++ b/test/Test_TIMap/5PO.mol2 @@ -0,0 +1,22 @@ +@MOLECULE +modXNA-5'-fragment:DPO:HEAD01:@HO5':HEAD01STRIP:NONE:TAIL01:@O5':TAIL01STRIP:@C5',H5'1,H5'2,H5'3 + 6 5 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 HO5' -0.4750 0.1600 -0.0250 HO 1 5PO 0.437230 + 2 O5' 0.9380 -0.2250 -0.7640 OH 1 5PO -0.640780 + 3 C5' 2.1230 -0.0970 -0.0430 CT 1 5PO -5.280519 + 4 H5'1 2.1370 -0.7580 0.8200 H1 1 5PO 1.651005 + 5 H5'2 2.2680 0.9190 0.3020 H1 1 5PO 1.651005 + 6 H5'3 2.9430 -0.3730 -0.7000 H1 1 5PO 1.651005 +@BOND + 1 1 2 1 + 2 2 3 1 + 3 3 4 1 + 4 3 5 1 + 5 3 6 1 +@SUBSTRUCTURE + 1 5PO 1 **** 0 **** **** diff --git a/test/Test_TIMap/A9Z.mol2 b/test/Test_TIMap/A9Z.mol2 new file mode 100644 index 0000000000..d509b36c21 --- /dev/null +++ b/test/Test_TIMap/A9Z.mol2 @@ -0,0 +1,78 @@ +@MOLECULE +modXNA-fragment:A9Z:HEAD01:@C5':HEAD01STRIP:@O1,H1:TAIL01:@O3':TAIL01STRIP:@H11:HEAD03:@C1':HEAD03STRIP:@C4,H6,H7,H8 + 33 34 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 O1 3.1370 -1.0420 -1.7190 OH 1 A9Z -0.679341 + 2 H1 3.7670 -1.2190 -2.4030 HO 1 A9Z 0.449849 + 3 C5' 3.7840 -0.9700 -0.4760 CI 1 A9Z 0.093465 + 4 H5' 4.6750 -0.3530 -0.5330 H1 1 A9Z 0.031050 + 5 H5'' 4.0770 -1.9640 -0.1400 H1 1 A9Z 0.031050 + 6 C4' 2.8410 -0.3550 0.5400 CT 1 A9Z 0.172623 + 7 H4' 3.2790 -0.4740 1.5280 H1 1 A9Z 0.018289 + 8 O4' 2.6910 1.0140 0.2570 OS 1 A9Z -0.489454 + 9 C1' 1.3330 1.4010 0.2740 CT 1 A9Z 0.413042 + 10 H1' 1.0290 1.6290 1.2870 H2 1 A9Z -0.007044 + 11 C4 1.1520 2.6120 -0.6170 CT 1 A9Z -0.414556 + 12 H6 0.1140 2.9250 -0.6200 HC 1 A9Z 0.143522 + 13 H7 1.4610 2.3870 -1.6320 HC 1 A9Z 0.143522 + 14 H8 1.7560 3.4340 -0.2500 HC 1 A9Z 0.143522 + 15 C3' 1.4230 -0.9620 0.5400 CT 1 A9Z 0.197117 + 16 H3' 1.3860 -1.9130 0.0220 H1 1 A9Z 0.035991 + 17 C2' 0.6340 0.1230 -0.1970 CT 1 A9Z -0.304034 + 18 H2' 0.8150 0.0060 -1.2570 H1 1 A9Z 0.152362 + 19 O3' 0.8980 -1.0800 1.8310 OS 1 A9Z -0.644079 + 20 H11 1.3250 -1.7870 2.2940 HO 1 A9Z 0.412277 + 21 N1' -0.7940 0.0880 0.0170 N* 1 A9Z 0.646711 + 22 N2' -1.3950 0.8720 0.8840 NB 1 A9Z -0.249810 + 23 N3' -2.6330 0.5970 0.8900 NB 1 A9Z -0.395140 + 24 C7' -2.8780 -0.4070 0.0130 CT 1 A9Z 0.349873 + 25 C8' -4.2570 -0.9510 -0.2130 CT 1 A9Z -0.134519 + 26 C6' -1.6900 -0.7390 -0.5550 CW 1 A9Z -0.515361 + 27 H6' -1.4160 -1.4630 -1.2920 H1 1 A9Z 0.278123 + 28 C9' -5.2010 0.0640 -0.8640 CT 1 A9Z -0.005132 + 29 H8' -4.1890 -1.8410 -0.8310 H1 1 A9Z 0.061954 + 30 H8'' -4.6690 -1.2650 0.7410 H1 1 A9Z 0.061954 + 31 H9'1 -6.1940 -0.3580 -0.9860 H1 1 A9Z 0.012460 + 32 H9'2 -5.2830 0.9550 -0.2530 H1 1 A9Z 0.012460 + 33 H9'3 -4.8370 0.3610 -1.8430 H1 1 A9Z 0.012460 +@BOND + 1 1 3 1 + 2 3 6 1 + 3 6 8 1 + 4 6 15 1 + 5 8 9 1 + 6 9 11 1 + 7 9 17 1 + 8 15 17 1 + 9 15 19 1 + 10 17 21 1 + 11 21 22 1 + 12 21 26 1 + 13 22 23 1 + 14 23 24 1 + 15 24 25 1 + 16 24 26 1 + 17 25 28 1 + 18 1 2 1 + 19 3 4 1 + 20 3 5 1 + 21 6 7 1 + 22 9 10 1 + 23 11 12 1 + 24 11 13 1 + 25 11 14 1 + 26 15 16 1 + 27 17 18 1 + 28 19 20 1 + 29 25 29 1 + 30 25 30 1 + 31 26 27 1 + 32 28 31 1 + 33 28 32 1 + 34 28 33 1 +@SUBSTRUCTURE + 1 A9Z 1 **** 0 **** **** diff --git a/test/Test_TIMap/DAA.mol2 b/test/Test_TIMap/DAA.mol2 new file mode 100644 index 0000000000..c1f4dc3ee4 --- /dev/null +++ b/test/Test_TIMap/DAA.mol2 @@ -0,0 +1,56 @@ +@MOLECULE +modXNA-fragment:DAA:HEAD01:@N9:HEAD01STRIP:@C1',H1',O4',HHO4,C2',H2'1,H2'3,H2'2 + 22 23 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 O4' -3.2440 -0.5840 0.8070 OH 1 DAA -0.745625 + 2 C1' -2.4070 0.3660 0.2200 CT 1 DAA 0.043100 + 3 H1' -2.2770 1.2210 0.8710 H2 1 DAA 0.183800 + 4 N9 -1.1080 -0.2670 0.1170 N* 1 DAA -0.026800 + 5 C8 -0.8260 -1.6030 -0.0450 C2 1 DAA 0.160700 + 6 H8 -1.6110 -2.3270 -0.0670 H5 1 DAA 0.187700 + 7 N7 0.4250 -1.8610 -0.1380 NB 1 DAA -0.617500 + 8 C5 1.0350 -0.6260 -0.0370 CB 1 DAA 0.072500 + 9 C6 2.3710 -0.2130 -0.0570 CA 1 DAA 0.689700 + 10 N6 3.3900 -1.0750 -0.2320 N2 1 DAA -0.912300 + 11 H61 4.3130 -0.7420 -0.0690 H 1 DAA 0.416700 + 12 H62 3.2170 -2.0510 -0.1560 H 1 DAA 0.416700 + 13 N1 2.6400 1.0800 0.0800 NC 1 DAA -0.762400 + 14 C2 1.6280 1.9270 0.2190 CQ 1 DAA 0.571600 + 15 H2 1.9070 2.9600 0.3220 H5 1 DAA 0.059800 + 16 N3 0.3410 1.6700 0.2480 NC 1 DAA -0.741700 + 17 C4 0.0970 0.3690 0.1220 CB 1 DAA 0.380000 + 18 C2' -2.9210 0.8200 -1.1380 CT 1 DAA -0.443298 + 19 H2'1 -3.0350 -0.0310 -1.8000 HC 1 DAA 0.124800 + 20 HHO4 -4.1330 -0.2580 0.8240 HO 1 DAA 0.478151 + 21 H2'2 -3.8850 1.3110 -1.0310 HC 1 DAA 0.124800 + 22 H2'3 -2.2310 1.5270 -1.5830 HC 1 DAA 0.124800 +@BOND + 1 1 2 1 + 2 2 4 1 + 3 2 18 1 + 4 4 5 1 + 5 4 17 1 + 6 5 7 1 + 7 7 8 1 + 8 8 9 1 + 9 8 17 1 + 10 9 10 1 + 11 9 13 1 + 12 13 14 1 + 13 14 16 1 + 14 16 17 1 + 15 1 20 1 + 16 2 3 1 + 17 5 6 1 + 18 10 11 1 + 19 10 12 1 + 20 14 15 1 + 21 18 19 1 + 22 18 21 1 + 23 18 22 1 +@SUBSTRUCTURE + 1 DAA 1 **** 0 **** **** diff --git a/test/Test_TIMap/DDD.mol2 b/test/Test_TIMap/DDD.mol2 new file mode 100644 index 0000000000..660dd989c3 --- /dev/null +++ b/test/Test_TIMap/DDD.mol2 @@ -0,0 +1,69 @@ +@MOLECULE +modXNA-fragment:DDD:HEAD01:@C5':HEAD01STRIP:@O5',P,OP1,OP2,OP3,H20:TAIL01:@O3':TAIL01STRIP:@P2,O12,O15,O16,H24:HEAD03:@C1':HEAD03STRIP:@C12,H21,H22,H23 + 29 29 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 P 3.3120 0.8260 0.2850 P 1 DDD 1.357695 + 2 OP1 3.6570 0.4210 1.6600 O2 1 DDD -0.844048 + 3 OP2 3.8210 2.0600 -0.3390 O2 1 DDD -0.844048 + 4 O5' 1.6640 0.7980 0.2240 OS 1 DDD -0.521429 + 5 O3' -2.2180 0.1650 -0.7670 OS 1 DDD -0.522241 + 6 C1' 0.2200 -2.4400 -0.1910 CT 1 DDD 0.383655 + 7 C2' -1.2130 -1.9410 -0.0230 CT 1 DDD -0.330047 + 8 C3' -1.0840 -0.4310 -0.2080 C7 1 DDD 0.254753 + 9 C4' 0.0840 -0.3150 -1.1950 CT 1 DDD 0.066591 + 10 C5' 0.9900 0.8930 -0.9940 CJ 1 DDD 0.148110 + 11 O4' 0.8360 -1.5150 -1.0720 OS 1 DDD -0.490259 + 12 C12 1.0120 -2.5340 1.1080 CT 1 DDD -0.293809 + 13 H1' 0.2430 -3.4090 -0.6840 H2 1 DDD -0.060446 + 14 H2' -1.8440 -2.3190 -0.8220 HC 1 DDD 0.116231 + 15 H2'' -1.6690 -2.2070 0.9200 HC 1 DDD 0.116231 + 16 H3' -0.8200 0.0460 0.7280 H1 1 DDD 0.052484 + 17 H4' -0.3370 -0.2720 -2.1970 H1 1 DDD 0.044462 + 18 H5' 1.6860 0.9680 -1.8200 H1 1 DDD 0.017730 + 19 H5'' 0.3690 1.7850 -1.0000 H1 1 DDD 0.017730 + 20 OP3 3.6120 -0.3960 -0.7630 OH 1 DDD -0.736711 + 21 H20 2.9190 -1.0440 -0.7120 HO 1 DDD 0.395040 + 22 H21 0.5560 -3.2640 1.7750 HC 1 DDD 0.087628 + 23 H22 1.0620 -1.5760 1.6110 HC 1 DDD 0.087628 + 24 H23 2.0320 -2.8460 0.9090 HC 1 DDD 0.087628 + 25 P2 -3.5100 0.5750 0.1540 P 1 DDD 1.280167 + 26 O12 -4.6060 0.8720 -0.7750 O2 1 DDD -0.822080 + 27 O15 -2.9840 1.9790 0.8140 OH 1 DDD -0.713044 + 28 O16 -3.6100 -0.3910 1.2730 O2 1 DDD -0.822080 + 29 H24 -2.6520 1.7830 1.6800 HO 1 DDD 0.397510 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 1 4 1 + 4 1 20 1 + 5 4 10 1 + 6 5 8 1 + 7 5 25 1 + 8 6 7 1 + 9 6 11 1 + 10 6 12 1 + 11 7 8 1 + 12 8 9 1 + 13 9 10 1 + 14 9 11 1 + 15 25 26 1 + 16 25 27 1 + 17 25 28 1 + 18 6 13 1 + 19 7 14 1 + 20 7 15 1 + 21 8 16 1 + 22 9 17 1 + 23 10 18 1 + 24 10 19 1 + 25 12 22 1 + 26 12 23 1 + 27 12 24 1 + 28 20 21 1 + 29 27 29 1 +@SUBSTRUCTURE + 1 DDD 1 **** 0 **** **** diff --git a/test/Test_TIMap/DPO.mol2 b/test/Test_TIMap/DPO.mol2 new file mode 100644 index 0000000000..fbde0e2584 --- /dev/null +++ b/test/Test_TIMap/DPO.mol2 @@ -0,0 +1,100 @@ +@MOLECULE +modXNA-fragment:DPO:HEAD01:@P:HEAD01STRIP:@O3',C3,H31,H32,H33:TAIL01:@O5':TAIL01STRIP:@C5',H5'1,H5'2,H5'3 + 13 12 1 0 0 +NUCLEIC_ACID +USER_CHARGES + + +@ATOM + 1 P -0.808000 -8.873000 -2.080000 P 1 DPO 1.1659 + 2 OP1 -2.029000 -8.583000 -1.290000 O2 1 DPO -0.7761 + 3 OP2 -0.196000 -10.198000 -1.860000 O2 1 DPO -0.7761 + 4 O5' 0.270000 -7.725000 -1.830000 OS 1 DPO -0.4954 + 5 C5' 1.364000 -7.578000 -2.770000 CJ 1 DPO 0.1500 + 6 H5'1 1.104000 -6.836000 -3.525000 H1 1 DPO 0.0182 + 7 H5'2 1.561000 -8.532000 -3.258000 H1 1 DPO 0.0182 + 8 O3' -1.195228 -8.803498 -3.527489 OS 1 DPO -0.5232 + 9 C3 -1.779497 -7.454439 -3.825261 CJ 1 DPO 0.1500 + 10 H5'3 2.257559 -7.253395 -2.236818 H1 1 DPO 0.0182 + 11 H31 -2.060879 -7.403911 -4.877056 H1 1 DPO 0.0181 + 12 H32 -1.040768 -6.682402 -3.609846 H1 1 DPO 0.0181 + 13 H33 -2.662006 -7.295315 -3.205566 H1 1 DPO 0.0181 +@BOND + 1 1 4 1 + 2 4 5 1 + 3 1 2 2 + 4 1 3 1 + 5 6 5 1 + 6 7 5 1 + 7 1 8 1 + 8 8 9 1 + 9 5 10 1 + 10 9 11 1 + 11 9 12 1 + 12 9 13 1 +@SUBSTRUCTURE +1 DPO 1 RESIDUE 1 S DA +@SET +GROUP:Nucleic_Acid_Groups:Sugar STATIC ATOMS +5 4 5 6 7 10 +GROUP:Nucleic_Acid_Groups:Phosphate STATIC ATOMS +4 1 2 3 8 +GROUP:Nucleic_Acid_Groups:Backbone STATIC ATOMS +9 1 2 3 4 5 6 7 8 10 +GROUP:Nucleic_Acid_Groups:Purine STATIC SUBSTRUCTURES +1 1 +GROUP:Nucleic_Acid_Groups:Base STATIC ATOMS +4 9 11 12 13 +GROUP:Nucleic_Acid_Groups:Pyrimidine STATIC SUBSTRUCTURES +1 1 +@ATOM_PROPERTY +Formal Charge +SciTegic.value.StringValue +1 0 +2 0 +3 -1 +4 0 +5 0 +6 0 +7 0 +8 0 +9 0 +10 0 +11 0 +12 0 +13 0 + +@ATOM_PROPERTY +pdb_tempFactor +SciTegic.value.DoubleValue +1 0.00 +2 0.00 +3 0.00 +4 0.00 +5 0.00 +6 0.00 +7 0.00 +8 0.00 +9 0.00 +10 0.00 +11 0.00 +12 0.00 +13 0.00 + +@ATOM_PROPERTY +pdb_occupancy +SciTegic.value.DoubleValue +1 1.00 +2 1.00 +3 1.00 +4 1.00 +5 1.00 +6 1.00 +7 1.00 +8 1.00 +9 1.00 +10 1.00 +11 1.00 +12 1.00 +13 1.00 + diff --git a/test/Test_TIMap/DTT.mol2 b/test/Test_TIMap/DTT.mol2 new file mode 100644 index 0000000000..a9aaa6d238 --- /dev/null +++ b/test/Test_TIMap/DTT.mol2 @@ -0,0 +1,55 @@ +@MOLECULE +modXNA-fragment:DTT:HEAD01:@N1:HEAD01STRIP:@C1',H1',O4',HHO4,C2',H2'1,H2'3,H2'2 + 22 22 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 C1' 2.1810 -0.1650 -0.1800 CT 1 TTT 0.067400 + 2 C2' 2.7740 -0.3090 1.2130 CT 1 TTT -0.415654 + 3 O4' 2.3730 -1.3160 -0.9490 OH 1 TTT -0.730006 + 4 O4 -3.1900 0.7190 0.0810 O 1 TTT -0.556300 + 5 C4 -2.0240 0.4570 0.0230 C 1 TTT 0.519400 + 6 N3 -1.0800 1.4670 -0.0700 NA 1 TTT -0.434000 + 7 C2 0.2820 1.3460 -0.1330 C 1 TTT 0.567700 + 8 O2 1.0060 2.3000 -0.1760 O 1 TTT -0.588100 + 9 N1 0.7370 0.0520 -0.1330 N* 1 TTT -0.023900 + 10 C6 -0.1350 -1.0170 -0.0340 CM 1 TTT -0.220900 + 11 C5 -1.4590 -0.8940 0.0440 CM 1 TTT 0.002500 + 12 C7 -2.4080 -2.0550 0.1420 CT 1 TTT -0.226900 + 13 H1' 2.5960 0.6980 -0.6740 H2 1 TTT 0.180400 + 14 H2'1 3.8480 -0.4630 1.1510 HC 1 TTT 0.116488 + 15 H2'2 2.5960 0.5910 1.7900 HC 1 TTT 0.116488 + 16 H3 -1.4260 2.4030 -0.0750 H 1 TTT 0.342000 + 17 H6 0.3440 -1.9740 -0.0420 H4 1 TTT 0.260700 + 18 H71 -3.0000 -1.9950 1.0480 HC 1 TTT 0.077000 + 19 H72 -3.1000 -2.0590 -0.6930 HC 1 TTT 0.077000 + 20 H73 -1.8680 -2.9950 0.1440 HC 1 TTT 0.077000 + 21 H2'3 2.3360 -1.1570 1.7300 HC 1 TTT 0.116488 + 22 HHO4 3.3010 -1.4750 -1.0520 HO 1 TTT 0.472667 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 1 9 1 + 4 4 5 1 + 5 5 6 1 + 6 5 11 1 + 7 6 7 1 + 8 7 8 1 + 9 7 9 1 + 10 9 10 1 + 11 10 11 1 + 12 11 12 1 + 13 1 13 1 + 14 2 14 1 + 15 2 15 1 + 16 2 21 1 + 17 3 22 1 + 18 6 16 1 + 19 10 17 1 + 20 12 18 1 + 21 12 19 1 + 22 12 20 1 +@SUBSTRUCTURE + 1 TTT 1 **** 0 **** **** diff --git a/test/Test_TIMap/ERN.lib b/test/Test_TIMap/ERN.lib new file mode 100644 index 0000000000..d8e56d108c --- /dev/null +++ b/test/Test_TIMap/ERN.lib @@ -0,0 +1,228 @@ +!!index array str + "ERN" +!entry.ERN.unit.atoms table str name str type int typex int resx int flags int seq int elmnt dbl chg + "N9" "N*" 0 1 131072 1 7 -0.026798 + "C8" "C2" 0 1 131072 2 6 0.160711 + "H8" "H5" 0 1 131072 3 1 0.187713 + "N7" "NB" 0 1 131072 4 7 -0.617458 + "C5" "CB" 0 1 131072 5 6 0.072505 + "C6" "CA" 0 1 131072 6 6 0.689747 + "N6" "N2" 0 1 131072 7 7 -0.912238 + "H61" "H" 0 1 131072 8 1 0.416728 + "H62" "H" 0 1 131072 9 1 0.416728 + "N1" "NC" 0 1 131072 10 7 -0.762348 + "C2" "CQ" 0 1 131072 11 6 0.571639 + "H2" "H5" 0 1 131072 12 1 0.059804 + "N3" "NC" 0 1 131072 13 7 -0.741649 + "C4" "CB" 0 1 131072 14 6 0.380026 + "P" "P" 0 1 131072 15 15 1.166744 + "OP1" "O2" 0 1 131072 16 8 -0.775638 + "OP2" "O2" 0 1 131072 17 8 -0.775638 + "O5'" "OS" 0 1 131072 18 8 -0.498667 + "C5'" "CJ" 0 1 131072 19 6 -0.006896 + "H5'" "H1" 0 1 131072 20 1 0.075442 + "H5''" "H1" 0 1 131072 21 1 0.075442 + "C4'" "CT" 0 1 131072 22 6 0.162991 + "H4'" "H1" 0 1 131072 23 1 0.117666 + "O4'" "OS" 0 1 131072 24 8 -0.368894 + "C1'" "CT" 0 1 131072 25 6 0.043124 + "C3'" "C7" 0 1 131072 26 6 0.071340 + "H3'" "H1" 0 1 131072 27 1 0.098555 + "C2'" "CT" 0 1 131072 28 6 -0.085352 + "H2'" "HC" 0 1 131072 29 1 0.071840 + "O3'" "OS" 0 1 131072 30 8 -0.522908 + "H1'" "H2" 0 1 131072 31 1 0.183902 + "H2''" "HC" 0 1 131072 32 1 0.071840 +!entry.ERN.unit.atomspertinfo table str pname str ptype int ptypex int pelmnt dbl pchg + "N9" "N*" 0 -1 0.0 + "C8" "C2" 0 -1 0.0 + "H8" "H5" 0 -1 0.0 + "N7" "NB" 0 -1 0.0 + "C5" "CB" 0 -1 0.0 + "C6" "CA" 0 -1 0.0 + "N6" "N2" 0 -1 0.0 + "H61" "H" 0 -1 0.0 + "H62" "H" 0 -1 0.0 + "N1" "NC" 0 -1 0.0 + "C2" "CQ" 0 -1 0.0 + "H2" "H5" 0 -1 0.0 + "N3" "NC" 0 -1 0.0 + "C4" "CB" 0 -1 0.0 + "P" "P" 0 -1 0.0 + "OP1" "O2" 0 -1 0.0 + "OP2" "O2" 0 -1 0.0 + "O5'" "OS" 0 -1 0.0 + "C5'" "CJ" 0 -1 0.0 + "H5'" "H1" 0 -1 0.0 + "H5''" "H1" 0 -1 0.0 + "C4'" "CT" 0 -1 0.0 + "H4'" "H1" 0 -1 0.0 + "O4'" "OS" 0 -1 0.0 + "C1'" "CT" 0 -1 0.0 + "C3'" "C7" 0 -1 0.0 + "H3'" "H1" 0 -1 0.0 + "C2'" "CT" 0 -1 0.0 + "H2'" "HC" 0 -1 0.0 + "O3'" "OS" 0 -1 0.0 + "H1'" "H2" 0 -1 0.0 + "H2''" "HC" 0 -1 0.0 +!entry.ERN.unit.boundbox array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.ERN.unit.childsequence single int + 2 +!entry.ERN.unit.connect array int + 15 + 30 +!entry.ERN.unit.connectivity table int atom1x int atom2x int flags + 1 2 1 + 1 14 1 + 1 25 1 + 2 4 1 + 2 3 1 + 4 5 1 + 5 6 1 + 5 14 1 + 6 7 1 + 6 10 1 + 7 8 1 + 7 9 1 + 10 11 1 + 11 13 1 + 11 12 1 + 13 14 1 + 15 18 1 + 15 16 1 + 15 17 1 + 18 19 1 + 19 22 1 + 19 20 1 + 19 21 1 + 22 24 1 + 22 26 1 + 22 23 1 + 24 25 1 + 25 28 1 + 25 31 1 + 26 28 1 + 26 30 1 + 26 27 1 + 28 29 1 + 28 32 1 +!entry.ERN.unit.hierarchy table str abovetype int abovex str belowtype int belowx + "U" 0 "R" 1 + "R" 1 "A" 1 + "R" 1 "A" 2 + "R" 1 "A" 3 + "R" 1 "A" 4 + "R" 1 "A" 5 + "R" 1 "A" 6 + "R" 1 "A" 7 + "R" 1 "A" 8 + "R" 1 "A" 9 + "R" 1 "A" 10 + "R" 1 "A" 11 + "R" 1 "A" 12 + "R" 1 "A" 13 + "R" 1 "A" 14 + "R" 1 "A" 15 + "R" 1 "A" 16 + "R" 1 "A" 17 + "R" 1 "A" 18 + "R" 1 "A" 19 + "R" 1 "A" 20 + "R" 1 "A" 21 + "R" 1 "A" 22 + "R" 1 "A" 23 + "R" 1 "A" 24 + "R" 1 "A" 25 + "R" 1 "A" 26 + "R" 1 "A" 27 + "R" 1 "A" 28 + "R" 1 "A" 29 + "R" 1 "A" 30 + "R" 1 "A" 31 + "R" 1 "A" 32 +!entry.ERN.unit.name single str + "Cpptraj" +!entry.ERN.unit.positions table dbl x dbl y dbl z + 2.456200 -0.350100 -0.692600 + 2.681900 0.897900 -1.227400 + 1.869200 1.536400 -1.560000 + 3.946900 1.238900 -1.292300 + 4.600500 0.119100 -0.743000 + 5.952100 -0.233600 -0.485700 + 7.007000 0.520000 -0.747300 + 7.921600 0.160900 -0.518900 + 6.878600 1.430300 -1.167300 + 6.272600 -1.405700 0.067200 + 5.286300 -2.241200 0.371400 + 5.581000 -3.183300 0.817500 + 3.979200 -2.061500 0.193500 + 3.700300 -0.847500 -0.371000 + -1.937600 3.244700 0.253100 + -2.098800 3.796900 -1.122900 + -1.079600 4.009200 1.203400 + -1.243500 1.790100 0.090200 + -1.901200 0.719300 -0.585100 + -2.046100 0.978200 -1.636000 + -2.883700 0.555200 -0.137300 + -1.114600 -0.602600 -0.509900 + -1.714100 -1.385100 -0.980300 + 0.120900 -0.474900 -1.212300 + 1.159600 -1.006600 -0.404800 + -0.787800 -1.027200 0.935900 + -1.336000 -0.422100 1.662200 + 0.710200 -0.757900 1.033600 + 0.882700 0.282800 1.313400 + -1.046000 -2.412200 1.148900 + 1.243600 -2.085100 -0.572000 + 1.204900 -1.431600 1.734800 +!entry.ERN.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x + 0 0 0 0 0 0 +!entry.ERN.unit.residues table str name int seq int childseq int startatomx str restype int imagingx + "ERN" 1 33 1 "n" 0 +!entry.ERN.unit.residuesPdbSequenceNumber array int + 0 +!entry.ERN.unit.solventcap array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.ERN.unit.velocities table dbl x dbl y dbl z + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 diff --git a/test/Test_TIMap/FLE.lib b/test/Test_TIMap/FLE.lib new file mode 100644 index 0000000000..6801aa2a56 --- /dev/null +++ b/test/Test_TIMap/FLE.lib @@ -0,0 +1,234 @@ +!!index array str + "FLE" +!entry.FLE.unit.atoms table str name str type int typex int resx int flags int seq int elmnt dbl chg + "N9" "N*" 0 1 131072 1 7 -0.025018 + "C8" "C5" 0 1 131072 2 6 0.201257 + "H8" "H5" 0 1 131072 3 1 0.155808 + "N7" "NB" 0 1 131072 4 7 -0.605312 + "C5" "CB" 0 1 131072 5 6 0.051669 + "C6" "CA" 0 1 131072 6 6 0.703195 + "N6" "N2" 0 1 131072 7 7 -0.898947 + "H61" "H" 0 1 131072 8 1 0.412847 + "H62" "H" 0 1 131072 9 1 0.412847 + "N1" "NC" 0 1 131072 10 7 -0.759007 + "C2" "CQ" 0 1 131072 11 6 0.589424 + "H2" "H5" 0 1 131072 12 1 0.047455 + "N3" "NC" 0 1 131072 13 7 -0.697409 + "C4" "CB" 0 1 131072 14 6 0.306300 + "P" "P" 0 1 131072 15 15 1.166744 + "OP1" "O2" 0 1 131072 16 8 -0.775638 + "OP2" "O2" 0 1 131072 17 8 -0.775638 + "O5'" "OS" 0 1 131072 18 8 -0.498667 + "C5'" "CI" 0 1 131072 19 6 0.055388 + "H5'" "H1" 0 1 131072 20 1 0.067399 + "H5''" "H1" 0 1 131072 21 1 0.067399 + "C4'" "CT" 0 1 131072 22 6 0.105714 + "H4'" "H1" 0 1 131072 23 1 0.116533 + "O4'" "OS" 0 1 131072 24 8 -0.357420 + "C1'" "CT" 0 1 131072 25 6 0.039109 + "C3'" "CT" 0 1 131072 26 6 0.200707 + "H3'" "H1" 0 1 131072 27 1 0.061046 + "C2'" "CT" 0 1 131072 28 6 0.066505 + "H2'" "H1" 0 1 131072 29 1 0.096482 + "O2'" "OH" 0 1 131072 30 8 -0.618434 + "O3'" "OS" 0 1 131072 31 8 -0.527064 + "H1'" "H2" 0 1 131072 32 1 0.199218 + "HO2'" "HO" 0 1 131072 33 1 0.415509 +!entry.FLE.unit.atomspertinfo table str pname str ptype int ptypex int pelmnt dbl pchg + "N9" "N*" 0 -1 0.0 + "C8" "C5" 0 -1 0.0 + "H8" "H5" 0 -1 0.0 + "N7" "NB" 0 -1 0.0 + "C5" "CB" 0 -1 0.0 + "C6" "CA" 0 -1 0.0 + "N6" "N2" 0 -1 0.0 + "H61" "H" 0 -1 0.0 + "H62" "H" 0 -1 0.0 + "N1" "NC" 0 -1 0.0 + "C2" "CQ" 0 -1 0.0 + "H2" "H5" 0 -1 0.0 + "N3" "NC" 0 -1 0.0 + "C4" "CB" 0 -1 0.0 + "P" "P" 0 -1 0.0 + "OP1" "O2" 0 -1 0.0 + "OP2" "O2" 0 -1 0.0 + "O5'" "OS" 0 -1 0.0 + "C5'" "CI" 0 -1 0.0 + "H5'" "H1" 0 -1 0.0 + "H5''" "H1" 0 -1 0.0 + "C4'" "CT" 0 -1 0.0 + "H4'" "H1" 0 -1 0.0 + "O4'" "OS" 0 -1 0.0 + "C1'" "CT" 0 -1 0.0 + "C3'" "CT" 0 -1 0.0 + "H3'" "H1" 0 -1 0.0 + "C2'" "CT" 0 -1 0.0 + "H2'" "H1" 0 -1 0.0 + "O2'" "OH" 0 -1 0.0 + "O3'" "OS" 0 -1 0.0 + "H1'" "H2" 0 -1 0.0 + "HO2'" "HO" 0 -1 0.0 +!entry.FLE.unit.boundbox array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.FLE.unit.childsequence single int + 2 +!entry.FLE.unit.connect array int + 15 + 31 +!entry.FLE.unit.connectivity table int atom1x int atom2x int flags + 1 2 1 + 1 14 1 + 1 25 1 + 2 4 1 + 2 3 1 + 4 5 1 + 5 6 1 + 5 14 1 + 6 7 1 + 6 10 1 + 7 8 1 + 7 9 1 + 10 11 1 + 11 13 1 + 11 12 1 + 13 14 1 + 15 18 1 + 15 16 1 + 15 17 1 + 18 19 1 + 19 22 1 + 19 20 1 + 19 21 1 + 22 24 1 + 22 26 1 + 22 23 1 + 24 25 1 + 25 28 1 + 25 32 1 + 26 28 1 + 26 31 1 + 26 27 1 + 28 30 1 + 28 29 1 + 30 33 1 +!entry.FLE.unit.hierarchy table str abovetype int abovex str belowtype int belowx + "U" 0 "R" 1 + "R" 1 "A" 1 + "R" 1 "A" 2 + "R" 1 "A" 3 + "R" 1 "A" 4 + "R" 1 "A" 5 + "R" 1 "A" 6 + "R" 1 "A" 7 + "R" 1 "A" 8 + "R" 1 "A" 9 + "R" 1 "A" 10 + "R" 1 "A" 11 + "R" 1 "A" 12 + "R" 1 "A" 13 + "R" 1 "A" 14 + "R" 1 "A" 15 + "R" 1 "A" 16 + "R" 1 "A" 17 + "R" 1 "A" 18 + "R" 1 "A" 19 + "R" 1 "A" 20 + "R" 1 "A" 21 + "R" 1 "A" 22 + "R" 1 "A" 23 + "R" 1 "A" 24 + "R" 1 "A" 25 + "R" 1 "A" 26 + "R" 1 "A" 27 + "R" 1 "A" 28 + "R" 1 "A" 29 + "R" 1 "A" 30 + "R" 1 "A" 31 + "R" 1 "A" 32 + "R" 1 "A" 33 +!entry.FLE.unit.name single str + "Cpptraj" +!entry.FLE.unit.positions table dbl x dbl y dbl z + -2.021500 -0.585900 -1.355200 + -1.515300 -1.823400 -1.686700 + -0.469400 -1.975600 -1.927700 + -2.398000 -2.792700 -1.681900 + -3.589700 -2.133000 -1.336700 + -4.942300 -2.529400 -1.155500 + -5.403300 -3.762000 -1.283500 + -6.387400 -3.924200 -1.134600 + -4.767800 -4.507800 -1.529800 + -5.889200 -1.644900 -0.835300 + -5.525900 -0.375500 -0.692800 + -6.315600 0.318700 -0.435800 + -4.308700 0.148700 -0.826900 + -3.370900 -0.791500 -1.147800 + 2.482500 -2.848700 -0.042000 + 3.207400 -3.002800 -1.335900 + 1.636800 -3.989900 0.410500 + 1.481400 -1.592600 -0.222500 + 1.981900 -0.275300 -0.420300 + 2.583800 -0.239900 -1.330600 + 2.618400 0.002900 0.422300 + 0.855800 0.764200 -0.540400 + 1.310500 1.752900 -0.624600 + 0.049200 0.551900 -1.698700 + -1.327900 0.731800 -1.362800 + -0.092800 0.788000 0.663800 + -0.310000 -0.232700 0.986100 + -1.352900 1.389500 0.027300 + -2.239600 1.100300 0.595900 + -1.291500 2.811900 -0.125200 + 0.387700 1.556500 1.768100 + -1.824700 1.385200 -2.085000 + -1.126600 3.214200 0.751700 +!entry.FLE.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x + 0 0 0 0 0 0 +!entry.FLE.unit.residues table str name int seq int childseq int startatomx str restype int imagingx + "FLE" 1 34 1 "n" 0 +!entry.FLE.unit.residuesPdbSequenceNumber array int + 0 +!entry.FLE.unit.solventcap array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.FLE.unit.velocities table dbl x dbl y dbl z + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 diff --git a/test/Test_TIMap/OME.mol2 b/test/Test_TIMap/OME.mol2 new file mode 100644 index 0000000000..81fcf40f28 --- /dev/null +++ b/test/Test_TIMap/OME.mol2 @@ -0,0 +1,61 @@ +@MOLECULE +modXNA-fragment:OME:HEAD01:@C5':HEAD01STRIP:@O1,H6:TAIL01:@O3':TAIL01STRIP:@H7:HEAD03:@C1':HEAD03STRIP:@C6,H9,H10,H11 + 25 25 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 O1 -2.7373 -0.2986 0.9909 OH 1 OME -0.658348 + 2 C5' -2.7524 0.3397 -0.2608 CI 1 OME 0.119115 + 3 H5' -3.0793 1.3567 -0.0940 H1 1 OME 0.042013 + 4 H5'' -3.4681 -0.1362 -0.9296 H1 1 OME 0.042013 + 5 C4' -1.3842 0.3812 -0.9232 CT 1 OME 0.116491 + 6 H4' -1.5101 0.7720 -1.9318 H1 1 OME 0.013915 + 7 O4' -0.5514 1.2469 -0.1927 OS 1 OME -0.474448 + 8 C1' 0.6960 0.6492 0.0801 CT 1 OME 0.283784 + 9 C3' -0.7976 -0.9215 -0.6469 CT 1 OME 0.127930 + 10 H3' -1.1941 -1.1239 0.3457 H1 1 OME 0.019321 + 11 C2' 0.6823 -0.7355 -0.6333 CT 1 OME 0.174570 + 12 H2' 1.4551 -1.4302 -0.3109 H1 1 OME 0.028124 + 13 O3' -1.1700 -2.0638 -1.0949 OS 1 OME -0.658945 + 14 H6 -2.9513 -1.2139 0.8920 HO 1 OME 0.409446 + 15 C6 1.2669 1.2326 1.3547 CT 1 OME -0.214186 + 16 H7 -1.4631 -2.8504 -0.6568 HO 1 OME 0.425460 + 17 H1' 1.3773 0.8180 -0.7490 H2 1 OME 0.018264 + 18 O2' 1.0039 -0.1767 -1.8583 OS 1 OME -0.414486 + 19 H9 2.2104 0.7535 1.5917 H1 1 OME 0.075428 + 20 H10 0.5774 1.0868 2.1784 H1 1 OME 0.075428 + 21 H11 1.4405 2.2973 1.2391 H1 1 OME 0.075428 + 22 CM2 0.1755 0.8586 -2.2875 CT 1 OME -0.017009 + 23 HM'1 0.4874 1.1245 -3.2888 H1 1 OME 0.066018 + 24 HM'2 0.2676 1.7327 -1.6505 H1 1 OME 0.066018 + 25 HM'3 -0.8669 0.5454 -2.3163 H1 1 OME 0.066018 +@BOND + 1 1 2 1 + 2 2 5 1 + 3 5 7 1 + 4 5 9 1 + 5 7 8 1 + 6 8 11 1 + 7 8 15 1 + 8 9 11 1 + 9 9 13 1 + 10 11 18 1 + 11 18 22 1 + 12 1 14 1 + 13 2 3 1 + 14 2 4 1 + 15 5 6 1 + 16 8 17 1 + 17 9 10 1 + 18 11 12 1 + 19 13 16 1 + 20 15 19 1 + 21 15 20 1 + 22 15 21 1 + 23 22 23 1 + 24 22 24 1 + 25 22 25 1 +@SUBSTRUCTURE + 1 OME 1 **** 0 **** **** diff --git a/test/Test_TIMap/PARENTS.txt b/test/Test_TIMap/PARENTS.txt new file mode 100644 index 0000000000..c4fec853e3 --- /dev/null +++ b/test/Test_TIMap/PARENTS.txt @@ -0,0 +1,30 @@ +# Example ModXNA parent fragments for the timap tests. +# +# These 11 fragments are consensus parents from the ModXNA library +# (https://github.com/modna-ff). They are test data, not files the timap +# command loads, and they are not installed under dat/. +# +# From this directory: +# parm DAA.mol2 name DAA +# timap analog template DAA out analog.lib +# +# Combined nucleotides are mapped to a user residue (for example FLE = rA +# onto ERN = dA), not to these fragments. Use naorder on that residue to +# freeze a full-nucleotide template. For amino acids, aaorder walks the +# residue in ff19SB amino19.lib order without reading that library. A +# smaller amino-acid library already in the test tree is +# test/Test_Build/amino12.lib. Small molecules: generate your own parent. +# +# file class notes +# ----- ----------------------- ------------------------------------------ +# 5PO backbone, no phosphate 5'-OH cap +# DPO backbone, phosphate canonical phosphodiester +# DAA purine base deoxyadenine nucleobase +# DTT pyrimidine base deoxythymine nucleobase +# TMC purine base + phosphate rare combined fragment; 1 corpus member +# UFB sugar, 2'-F +# A9Z sugar, 2'-HH bulky analog; prefer a real d-ribose if you have one +# DDD sugar, 2'-HH + phosphate rare combined fragment; 1 corpus member +# RUR sugar, 2'-OH ribose +# OME sugar, 2'-OMe +# RBO sugar, 2'-OR 2'-O-benzyl, not plain ribose diff --git a/test/Test_TIMap/RBO.mol2 b/test/Test_TIMap/RBO.mol2 new file mode 100644 index 0000000000..bc68ca73d4 --- /dev/null +++ b/test/Test_TIMap/RBO.mol2 @@ -0,0 +1,82 @@ +@MOLECULE +modXNA-fragment:RBO:HEAD01:@C5':HEAD01STRIP:@O1,H6:TAIL01:@O3':TAIL01STRIP:@H7:HEAD03:@C1':HEAD03STRIP:@C6,H9,H10,H11 + 35 36 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 O1 -4.7030 -0.8700 0.0080 OH 1 RBO -0.686540 + 2 C5' -3.6320 -1.6730 -0.3980 CI 1 RBO 0.218167 + 3 H5' -3.8010 -2.0640 -1.3980 H1 1 RBO 0.013241 + 4 H5'' -3.5780 -2.5090 0.2870 H1 1 RBO 0.013241 + 5 C4' -2.3190 -0.9140 -0.3690 CT 1 RBO 0.108579 + 6 H4' -1.5290 -1.5420 -0.7760 H1 1 RBO -0.000075 + 7 O4' -2.4690 0.2550 -1.1450 OS 1 RBO -0.480597 + 8 C1' -1.6410 1.3050 -0.6710 CT 1 RBO 0.209892 + 9 C3' -1.9000 -0.3920 1.0020 CT 1 RBO 0.221414 + 10 H3' -2.7760 -0.0250 1.5190 H1 1 RBO 0.062201 + 11 C2' -0.9930 0.7790 0.6200 CT 1 RBO -0.020677 + 12 H2' -0.9330 1.5350 1.3970 H1 1 RBO 0.075189 + 13 O2' 0.2820 0.2340 0.3890 OS 1 RBO -0.261827 + 14 O3' -1.2990 -1.3380 1.8210 OS 1 RBO -0.675377 + 15 H6 -4.7270 -0.1150 -0.5650 HO 1 RBO 0.426097 + 16 C6 -2.4890 2.5530 -0.4810 CT 1 RBO -0.314105 + 17 H7 -0.3920 -1.4250 1.5500 HO 1 RBO 0.395637 + 18 H1' -0.8730 1.4940 -1.4120 H2 1 RBO 0.067178 + 19 H9 -1.8710 3.4040 -0.2110 HC 1 RBO 0.110105 + 20 H10 -3.2300 2.4060 0.2990 HC 1 RBO 0.110105 + 21 H11 -3.0110 2.7890 -1.4010 HC 1 RBO 0.110105 + 22 C6' 1.3540 1.1380 0.4220 CT 1 RBO 0.013309 + 23 H6' 1.1860 1.9370 -0.2980 H1 1 RBO 0.058543 + 24 H6'' 1.4190 1.5960 1.4050 H1 1 RBO 0.058543 + 25 C7' 2.6380 0.4150 0.0980 CA 1 RBO 0.076537 + 26 C82' 2.7110 -0.4330 -1.0030 CA 1 RBO -0.150705 + 27 H82' 1.8350 -0.5980 -1.6040 HA 1 RBO 0.140745 + 28 C92' 3.8960 -1.0710 -1.3190 CA 1 RBO -0.087769 + 29 H92' 3.9390 -1.7270 -2.1700 HA 1 RBO 0.135949 + 30 C10' 5.0280 -0.8660 -0.5420 CA 1 RBO -0.118207 + 31 H10' 5.9500 -1.3630 -0.7900 HA 1 RBO 0.139145 + 32 C91' 4.9640 -0.0240 0.5520 CA 1 RBO -0.087769 + 33 H91' 5.8350 0.1360 1.1630 HA 1 RBO 0.135949 + 34 C81' 3.7720 0.6100 0.8720 CA 1 RBO -0.150705 + 35 H81' 3.7290 1.2560 1.7320 HA 1 RBO 0.140745 +@BOND + 1 1 2 1 + 2 2 5 1 + 3 5 7 1 + 4 5 9 1 + 5 7 8 1 + 6 8 11 1 + 7 8 16 1 + 8 9 11 1 + 9 9 14 1 + 10 11 13 1 + 11 13 22 1 + 12 22 25 1 + 13 25 26 1 + 14 25 34 1 + 15 26 28 1 + 16 28 30 1 + 17 30 32 1 + 18 32 34 1 + 19 1 15 1 + 20 2 3 1 + 21 2 4 1 + 22 5 6 1 + 23 8 18 1 + 24 9 10 1 + 25 11 12 1 + 26 14 17 1 + 27 16 19 1 + 28 16 20 1 + 29 16 21 1 + 30 22 23 1 + 31 22 24 1 + 32 26 27 1 + 33 28 29 1 + 34 30 31 1 + 35 32 33 1 + 36 34 35 1 +@SUBSTRUCTURE + 1 RBO 1 **** 0 **** **** diff --git a/test/Test_TIMap/RUR.mol2 b/test/Test_TIMap/RUR.mol2 new file mode 100644 index 0000000000..d14f7d7fa2 --- /dev/null +++ b/test/Test_TIMap/RUR.mol2 @@ -0,0 +1,58 @@ +@MOLECULE +modXNA-fragment:RUR:HEAD01:@C5':HEAD01STRIP:@O1,H6:TAIL01:@O3':TAIL01STRIP:@H7:HEAD03:@C1':HEAD03STRIP:@C6,H9,H10,H11 + 24 23 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 O1 -2.3050 -1.6420 -0.0240 OH 1 RUR -0.710314 + 2 C5' -1.6370 -0.8720 0.9380 CI 1 RUR 0.222004 + 3 H5' -0.9510 -1.4870 1.5170 H1 1 RUR 0.010411 + 4 H5'' -2.3890 -0.4930 1.6150 H1 1 RUR 0.010411 + 5 C4' -0.8740 0.2810 0.3020 CT 1 RUR 0.092794 + 6 H4' -0.4800 0.9130 1.0940 H1 1 RUR 0.025391 + 7 O4' 0.1700 -0.2280 -0.4980 OS 1 RUR -0.642531 + 8 C1' 1.4480 -0.3580 0.0820 CT 1 RUR 0.426540 + 9 C3' -1.7480 1.1200 -0.6280 CT 1 RUR 0.222004 + 10 H3' -1.9270 0.5600 -1.5390 H1 1 RUR 0.007912 + 11 C2' 2.2210 0.9420 -0.1330 CT 1 RUR 0.007159 + 12 H2' 2.3890 1.0840 -1.1920 H1 1 RUR 0.096729 + 13 O2' 3.4820 0.9180 0.4790 OH 1 RUR -0.653987 + 14 O3' -2.9470 1.5070 -0.0210 OS 1 RUR -0.710314 + 15 H6 -1.6680 -1.9440 -0.6590 HO 1 RUR 0.455382 + 16 C6 2.1400 -1.5500 -0.5560 CT 1 RUR -0.282455 + 17 H7 -3.5620 0.7870 -0.0760 HO 1 RUR 0.455382 + 18 H1' 1.3500 -0.5340 1.1530 H2 1 RUR -0.065373 + 19 H9 3.1530 -1.6450 -0.1860 HC 1 RUR 0.099165 + 20 H10 2.1750 -1.4250 -1.6330 HC 1 RUR 0.099165 + 21 H11 1.5990 -2.4650 -0.3380 HC 1 RUR 0.099165 + 22 H3'' -1.2160 2.0230 -0.8980 H1 1 RUR 0.010411 + 23 H2'' 1.6330 1.7870 0.2200 H1 1 RUR 0.096729 + 24 HO2' 3.3870 1.0240 1.4150 HO 1 RUR 0.408443 +@BOND + 1 1 2 1 + 2 2 5 1 + 3 5 7 1 + 4 5 9 1 + 5 7 8 1 + 6 8 11 1 + 7 8 16 1 + 8 9 14 1 + 9 11 13 1 + 10 1 15 1 + 11 2 3 1 + 12 2 4 1 + 13 5 6 1 + 14 8 18 1 + 15 9 10 1 + 16 9 22 1 + 17 11 12 1 + 18 11 23 1 + 19 13 24 1 + 20 14 17 1 + 21 16 19 1 + 22 16 20 1 + 23 16 21 1 +@SUBSTRUCTURE + 1 RUR 1 **** 0 **** **** diff --git a/test/Test_TIMap/RunTest.sh b/test/Test_TIMap/RunTest.sh new file mode 100755 index 0000000000..ae69922638 --- /dev/null +++ b/test/Test_TIMap/RunTest.sh @@ -0,0 +1,229 @@ +#!/bin/bash + +. ../MasterTest.sh + +CleanFiles timap.in \ + fle_to_ern.map ern_naorder.map \ + phenol_to_benzene.map sec_to_cys.map \ + daa_identity.map \ + cys_aaorder.map cys_shuffled.aaorder.map \ + fle.sorted.lib FLE.sorted.lib \ + series_parent.dat \ + series_map_benzene.map series_map_phenol.map \ + mix_pdb_mol2.map path_mix.map \ + phenol_aligned.mol2 phenol_aligned.pdb phenol_mol2_out.map \ + phenol_mol2_rt.map phenol_pdb_rt.map \ + fle_ern.0.mol2 fle_ern.1.mol2 fle_ern.0.lib fle_ern.1.lib \ + fle_ern.scmask fle_ern.atoms fle_ern.map \ + phenol_bnz.0.mol2 phenol_bnz.1.mol2 phenol_bnz.0.lib phenol_bnz.1.lib \ + phenol_bnz.scmask phenol_bnz.atoms phenol_bnz.map \ + sec_cys.0.mol2 sec_cys.1.mol2 sec_cys.0.lib sec_cys.1.lib \ + sec_cys.scmask sec_cys.atoms sec_cys.map + +TESTNAME='TI map tests' +Requires maxthreads 1 + +INPUT='-i timap.in' + +# Nucleic acids: unmodified rA (FLE) onto unmodified dA (ERN). +# Shared atoms keep ERN indices; 2'-OH is an insertion; ERN H2'' is unmatched. +UNITNAME='Nucleic acid: FLE (rA) onto ERN (dA)' +CheckFor maxthreads 1 +if [ $? -eq 0 ] ; then + cat > timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in < timap.in <MOLECULE +modXNA-fragment:TMC:HEAD01:@N3:HEAD01STRIP:@C1,C2,O1,H1,H2,H3,H7,H8 + 40 41 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 C1 -6.0900 -0.6500 0.1500 CT 1 TMC 0.299384 + 2 C2 -6.1670 -2.0190 0.8110 CT 1 TMC -0.338476 + 3 O1 -7.0100 0.2470 0.6940 OH 1 TMC -0.734852 + 4 N1 -0.3990 1.1790 -0.0370 N2 1 TMC -0.645719 + 5 C3 -1.4390 0.4620 -0.4960 CA 1 TMC 0.584071 + 6 C4 -2.6930 0.5090 0.1320 CB 1 TMC -0.082278 + 7 N2 -3.1610 1.1710 1.2490 NB 1 TMC -0.621722 + 8 C5 -4.3910 0.8200 1.3310 CK 1 TMC 0.204850 + 9 N3 -4.7900 -0.0380 0.3350 N* 1 TMC -0.139760 + 10 C6 -3.6890 -0.2440 -0.4410 CB 1 TMC 0.396792 + 11 N4 -3.5750 -1.0120 -1.5230 NC 1 TMC -0.820978 + 12 C7 -2.3600 -0.9750 -2.0150 CQ 1 TMC 0.515957 + 13 N5 -1.3030 -0.3000 -1.5750 NC 1 TMC -0.847307 + 14 H1 -6.2310 -0.7480 -0.9190 H2 1 TMC 0.068649 + 15 H2 -7.1420 -2.4660 0.6390 HC 1 TMC 0.126559 + 16 H3 -5.4160 -2.6790 0.3930 HC 1 TMC 0.126559 + 17 C8 0.9590 1.0530 -0.5280 CT 1 TMC 0.128351 + 18 H4 -0.5430 1.6270 0.8410 H 1 TMC 0.344693 + 19 H5 -5.0820 1.1520 2.0750 H5 1 TMC 0.155797 + 20 H6 -2.1860 -1.5720 -2.8930 H5 1 TMC 0.040599 + 21 H7 -6.0080 -1.9280 1.8790 HC 1 TMC 0.126559 + 22 H8 -7.8760 -0.1360 0.6760 HO 1 TMC 0.476626 + 23 C9 1.7530 2.3120 -0.1930 CT 1 TMC 0.025710 + 24 C10 3.2020 2.2650 -0.6980 CT 1 TMC -0.123482 + 25 C11 4.0540 1.2550 0.0210 CM 1 TMC 0.019055 + 26 C12 4.7310 0.2790 -0.5610 CM 1 TMC -0.351725 + 27 P1 5.7790 -0.8650 0.3090 P 1 TMC 1.207689 + 28 H9 1.4390 0.1730 -0.1080 H1 1 TMC 0.053905 + 29 H10 0.9150 0.9110 -1.5970 H1 1 TMC 0.053905 + 30 H11 1.2530 3.1670 -0.6360 HC 1 TMC 0.014755 + 31 H12 1.7540 2.4740 0.8840 HC 1 TMC 0.014755 + 32 H13 3.2190 2.0740 -1.7660 HC 1 TMC 0.078356 + 33 H14 3.6480 3.2460 -0.5470 HC 1 TMC 0.078356 + 34 H15 4.1010 1.3590 1.0920 HA 1 TMC 0.124936 + 35 H16 4.7000 0.1450 -1.6290 HA 1 TMC 0.209140 + 36 O2 7.2080 -0.8580 0.0140 O2 1 TMC -0.745756 + 37 O3 5.0890 -2.2670 -0.0050 OH 1 TMC -0.676900 + 38 O4 5.4800 -0.5870 1.8530 OH 1 TMC -0.676900 + 39 H17 6.2900 -0.5120 2.3470 HO 1 TMC 0.463101 + 40 H18 5.7230 -2.9760 -0.0410 HO 1 TMC 0.463101 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 1 9 1 + 4 4 5 1 + 5 4 17 1 + 6 5 6 1 + 7 5 13 1 + 8 6 7 1 + 9 6 10 1 + 10 7 8 1 + 11 8 9 1 + 12 9 10 1 + 13 10 11 1 + 14 11 12 1 + 15 12 13 1 + 16 17 23 1 + 17 23 24 1 + 18 24 25 1 + 19 25 26 1 + 20 26 27 1 + 21 27 36 1 + 22 27 37 1 + 23 27 38 1 + 24 1 14 1 + 25 2 15 1 + 26 2 16 1 + 27 2 21 1 + 28 3 22 1 + 29 4 18 1 + 30 8 19 1 + 31 12 20 1 + 32 17 28 1 + 33 17 29 1 + 34 23 30 1 + 35 23 31 1 + 36 24 32 1 + 37 24 33 1 + 38 25 34 1 + 39 26 35 1 + 40 37 40 1 + 41 38 39 1 +@SUBSTRUCTURE + 1 TMC 1 **** 0 **** **** diff --git a/test/Test_TIMap/UFB.mol2 b/test/Test_TIMap/UFB.mol2 new file mode 100644 index 0000000000..d887dd4c9a --- /dev/null +++ b/test/Test_TIMap/UFB.mol2 @@ -0,0 +1,61 @@ +@MOLECULE +modXNA-fragment:UFB:HEAD01:@C5':HEAD01STRIP:@O3,H10:TAIL01:@O3':TAIL01STRIP:@H4:HEAD03:@C1':HEAD03STRIP:@C2,H11,H12,H13 + 25 25 1 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 C1' -1.4540 0.1120 0.7800 CT 1 UFB 0.276778 + 2 C2 -2.6610 0.9730 1.0880 CT 1 UFB -0.434204 + 3 C2' -1.3440 -0.3890 -0.6510 CT 1 UFB 0.102100 + 4 F2' -2.1090 -1.5080 -0.8650 F 1 UFB -0.271425 + 5 C3' 0.1420 -0.7060 -0.7750 CT 1 UFB 0.076654 + 6 O3' 0.3880 -1.9970 -0.2850 OS 1 UFB -0.646388 + 7 C4' 0.7910 0.3400 0.1720 CT 1 UFB 0.206737 + 8 O6' 1.4320 1.3570 -0.5070 OS 1 UFB -0.308597 + 9 C6' 0.6860 2.2240 -1.3120 CT 1 UFB 0.101541 + 10 C5' 1.8230 -0.2500 1.1320 CI 1 UFB -0.007930 + 11 O3 2.8540 -0.9240 0.4810 OH 1 UFB -0.647820 + 12 O4' -0.2760 0.8680 0.9290 OS 1 UFB -0.467153 + 13 H1' -1.4340 -0.7500 1.4420 H2 1 UFB 0.037811 + 14 H2' -1.6610 0.3630 -1.3590 H1 1 UFB 0.078816 + 15 H3' 0.5260 -0.6080 -1.7820 H1 1 UFB 0.079384 + 16 H4 -0.2490 -2.5950 -0.6540 HO 1 UFB 0.445843 + 17 H6'1 1.3800 2.9630 -1.6860 HC 1 UFB 0.022519 + 18 H6'2 -0.0920 2.7220 -0.7470 HC 1 UFB 0.022519 + 19 H6'3 0.2430 1.7140 -2.1640 HC 1 UFB 0.022519 + 20 H5' 2.2600 0.5800 1.6670 H1 1 UFB 0.129040 + 21 H5'' 1.3120 -0.8880 1.8420 H1 1 UFB 0.129040 + 22 H10 2.5100 -1.7490 0.1650 HO 1 UFB 0.438283 + 23 H11 -2.6250 1.3180 2.1140 HC 1 UFB 0.162427 + 24 H12 -3.5730 0.4030 0.9480 HC 1 UFB 0.162427 + 25 H13 -2.6880 1.8410 0.4370 HC 1 UFB 0.162427 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 1 12 1 + 4 3 4 1 + 5 3 5 1 + 6 5 6 1 + 7 5 7 1 + 8 7 8 1 + 9 7 10 1 + 10 7 12 1 + 11 8 9 1 + 12 10 11 1 + 13 1 13 1 + 14 2 23 1 + 15 2 24 1 + 16 2 25 1 + 17 3 14 1 + 18 5 15 1 + 19 6 16 1 + 20 9 17 1 + 21 9 18 1 + 22 9 19 1 + 23 10 20 1 + 24 10 21 1 + 25 11 22 1 +@SUBSTRUCTURE + 1 UFB 1 **** 0 **** **** diff --git a/test/Test_TIMap/benzene.mol2 b/test/Test_TIMap/benzene.mol2 new file mode 100644 index 0000000000..7688f484c8 --- /dev/null +++ b/test/Test_TIMap/benzene.mol2 @@ -0,0 +1,35 @@ +@MOLECULE +benzene TI parent + 12 12 1 0 0 +SMALL +NO_CHARGES + + +@ATOM + 1 C1 1.2000 0.0000 0.0000 C 1 BNZ 0.000000 + 2 H1 2.2800 0.0000 0.0000 H 1 BNZ 0.000000 + 3 C2 0.6000 1.0390 0.0000 C 1 BNZ 0.000000 + 4 H2 1.1400 1.9740 0.0000 H 1 BNZ 0.000000 + 5 C3 -0.6000 1.0390 0.0000 C 1 BNZ 0.000000 + 6 H3 -1.1400 1.9740 0.0000 H 1 BNZ 0.000000 + 7 C4 -1.2000 0.0000 0.0000 C 1 BNZ 0.000000 + 8 H4 -2.2800 0.0000 0.0000 H 1 BNZ 0.000000 + 9 C5 -0.6000 -1.0390 0.0000 C 1 BNZ 0.000000 + 10 H5 -1.1400 -1.9740 0.0000 H 1 BNZ 0.000000 + 11 C6 0.6000 -1.0390 0.0000 C 1 BNZ 0.000000 + 12 H6 1.1400 -1.9740 0.0000 H 1 BNZ 0.000000 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 3 4 1 + 4 3 5 1 + 5 5 6 1 + 6 5 7 1 + 7 7 8 1 + 8 7 9 1 + 9 9 10 1 + 10 9 11 1 + 11 11 12 1 + 12 11 1 1 +@SUBSTRUCTURE + 1 BNZ 1 **** 0 **** **** diff --git a/test/Test_TIMap/benzene.pdb b/test/Test_TIMap/benzene.pdb new file mode 100644 index 0000000000..6473fc9910 --- /dev/null +++ b/test/Test_TIMap/benzene.pdb @@ -0,0 +1,25 @@ +ATOM 1 C1 BNZ A 1 1.200 0.000 0.000 1.00 0.00 C +ATOM 2 H1 BNZ A 1 2.280 0.000 0.000 1.00 0.00 H +ATOM 3 C2 BNZ A 1 0.600 1.039 0.000 1.00 0.00 C +ATOM 4 H2 BNZ A 1 1.140 1.974 0.000 1.00 0.00 H +ATOM 5 C3 BNZ A 1 -0.600 1.039 0.000 1.00 0.00 C +ATOM 6 H3 BNZ A 1 -1.140 1.974 0.000 1.00 0.00 H +ATOM 7 C4 BNZ A 1 -1.200 0.000 0.000 1.00 0.00 C +ATOM 8 H4 BNZ A 1 -2.280 0.000 0.000 1.00 0.00 H +ATOM 9 C5 BNZ A 1 -0.600 -1.039 0.000 1.00 0.00 C +ATOM 10 H5 BNZ A 1 -1.140 -1.974 0.000 1.00 0.00 H +ATOM 11 C6 BNZ A 1 0.600 -1.039 0.000 1.00 0.00 C +ATOM 12 H6 BNZ A 1 1.140 -1.974 0.000 1.00 0.00 H +CONECT 1 2 +CONECT 1 3 +CONECT 3 4 +CONECT 3 5 +CONECT 5 6 +CONECT 5 7 +CONECT 7 8 +CONECT 7 9 +CONECT 9 10 +CONECT 9 11 +CONECT 11 12 +CONECT 11 1 +END diff --git a/test/Test_TIMap/cys.lib b/test/Test_TIMap/cys.lib new file mode 100644 index 0000000000..a20e788365 --- /dev/null +++ b/test/Test_TIMap/cys.lib @@ -0,0 +1,99 @@ +!!index array str + "CYS" +!entry.CYS.unit.atoms table str name str type int typex int resx int flags int seq int elmnt dbl chg + "N" "N" 0 1 131072 1 7 0.000000 + "H" "H" 0 1 131072 2 1 0.000000 + "CA" "C" 0 1 131072 3 6 0.000000 + "HA" "H" 0 1 131072 4 1 0.000000 + "CB" "C" 0 1 131072 5 6 0.000000 + "HB2" "H" 0 1 131072 6 1 0.000000 + "HB3" "H" 0 1 131072 7 1 0.000000 + "SG" "S" 0 1 131072 8 16 0.000000 + "HG" "H" 0 1 131072 9 1 0.000000 + "C" "C" 0 1 131072 10 6 0.000000 + "O" "O" 0 1 131072 11 8 0.000000 +!entry.CYS.unit.atomspertinfo table str pname str ptype int ptypex int pelmnt dbl pchg + "N" "N" 0 -1 0.0 + "H" "H" 0 -1 0.0 + "CA" "C" 0 -1 0.0 + "HA" "H" 0 -1 0.0 + "CB" "C" 0 -1 0.0 + "HB2" "H" 0 -1 0.0 + "HB3" "H" 0 -1 0.0 + "SG" "S" 0 -1 0.0 + "HG" "H" 0 -1 0.0 + "C" "C" 0 -1 0.0 + "O" "O" 0 -1 0.0 +!entry.CYS.unit.boundbox array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.CYS.unit.childsequence single int + 2 +!entry.CYS.unit.connect array int + 0 + 0 +!entry.CYS.unit.connectivity table int atom1x int atom2x int flags + 1 2 1 + 1 3 1 + 3 4 1 + 3 5 1 + 3 10 1 + 5 6 1 + 5 7 1 + 5 8 1 + 8 9 1 + 10 11 1 +!entry.CYS.unit.hierarchy table str abovetype int abovex str belowtype int belowx + "U" 0 "R" 1 + "R" 1 "A" 1 + "R" 1 "A" 2 + "R" 1 "A" 3 + "R" 1 "A" 4 + "R" 1 "A" 5 + "R" 1 "A" 6 + "R" 1 "A" 7 + "R" 1 "A" 8 + "R" 1 "A" 9 + "R" 1 "A" 10 + "R" 1 "A" 11 +!entry.CYS.unit.name single str + "CYS" +!entry.CYS.unit.positions table dbl x dbl y dbl z + 0.000000 0.000000 0.000000 + 0.000000 0.000000 1.010000 + 1.450000 0.000000 -0.500000 + 1.450000 0.890000 -1.090000 + 2.200000 -1.200000 0.200000 + 1.700000 -2.100000 -0.100000 + 3.200000 -1.100000 -0.200000 + 2.400000 -1.400000 2.000000 + 3.500000 -1.100000 2.300000 + 2.200000 1.200000 0.200000 + 3.400000 1.200000 0.400000 +!entry.CYS.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x + 0 0 0 0 0 0 +!entry.CYS.unit.residues table str name int seq int childseq int startatomx str restype int imagingx + "CYS" 1 12 1 "?" 0 +!entry.CYS.unit.residuesPdbSequenceNumber array int + 0 +!entry.CYS.unit.solventcap array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.CYS.unit.velocities table dbl x dbl y dbl z + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 diff --git a/test/Test_TIMap/cys.mol2 b/test/Test_TIMap/cys.mol2 new file mode 100644 index 0000000000..db8c048a79 --- /dev/null +++ b/test/Test_TIMap/cys.mol2 @@ -0,0 +1,32 @@ +@MOLECULE +cysteine TI parent + 11 10 1 0 0 +SMALL +NO_CHARGES + + +@ATOM + 1 N 0.0000 0.0000 0.0000 N 1 CYS 0.000000 + 2 H 0.0000 0.0000 1.0100 H 1 CYS 0.000000 + 3 CA 1.4500 0.0000 -0.5000 C 1 CYS 0.000000 + 4 HA 1.4500 0.8900 -1.0900 H 1 CYS 0.000000 + 5 CB 2.2000 -1.2000 0.2000 C 1 CYS 0.000000 + 6 HB2 1.7000 -2.1000 -0.1000 H 1 CYS 0.000000 + 7 HB3 3.2000 -1.1000 -0.2000 H 1 CYS 0.000000 + 8 SG 2.4000 -1.4000 2.0000 S 1 CYS 0.000000 + 9 HG 3.5000 -1.1000 2.3000 H 1 CYS 0.000000 + 10 C 2.2000 1.2000 0.2000 C 1 CYS 0.000000 + 11 O 3.4000 1.2000 0.4000 O 1 CYS 0.000000 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 3 4 1 + 4 3 5 1 + 5 3 10 1 + 6 5 6 1 + 7 5 7 1 + 8 5 8 1 + 9 8 9 1 + 10 10 11 1 +@SUBSTRUCTURE + 1 CYS 1 **** 0 **** **** diff --git a/test/Test_TIMap/cys_aaorder.map.save b/test/Test_TIMap/cys_aaorder.map.save new file mode 100644 index 0000000000..d4f1b0e1f5 --- /dev/null +++ b/test/Test_TIMap/cys_aaorder.map.save @@ -0,0 +1,15 @@ +# timap tgt='CYS[CYS]' template='CYS[CYS]' +# kind tgt=amino template=amino +# mapped= 11 insertion= 0 unmapped_template= 0 n_tgt= 11 n_tpl= 11 +#Out TgtAt TgtName TplAt TplName Ins + 1 1 N 1 N 0 + 2 2 H 2 H 0 + 3 3 CA 3 CA 0 + 4 4 HA 4 HA 0 + 5 5 CB 5 CB 0 + 6 6 HB2 6 HB2 0 + 7 7 HB3 7 HB3 0 + 8 8 SG 8 SG 0 + 9 9 HG 9 HG 0 + 10 10 C 10 C 0 + 11 11 O 11 O 0 diff --git a/test/Test_TIMap/cys_shuffled.aaorder.map.save b/test/Test_TIMap/cys_shuffled.aaorder.map.save new file mode 100644 index 0000000000..0a7555b7c4 --- /dev/null +++ b/test/Test_TIMap/cys_shuffled.aaorder.map.save @@ -0,0 +1,15 @@ +# timap tgt='CYS[CYS]' template='CYS[CYS]' +# kind tgt=amino template=amino +# mapped= 11 insertion= 0 unmapped_template= 0 n_tgt= 11 n_tpl= 11 +#Out TgtAt TgtName TplAt TplName Ins + 1 3 N 3 N 0 + 2 4 H 4 H 0 + 3 5 CA 5 CA 0 + 4 6 HA 6 HA 0 + 5 7 CB 7 CB 0 + 6 8 HB2 8 HB2 0 + 7 9 HB3 9 HB3 0 + 8 10 SG 10 SG 0 + 9 11 HG 11 HG 0 + 10 1 C 1 C 0 + 11 2 O 2 O 0 diff --git a/test/Test_TIMap/cys_shuffled.lib b/test/Test_TIMap/cys_shuffled.lib new file mode 100644 index 0000000000..d360c46806 --- /dev/null +++ b/test/Test_TIMap/cys_shuffled.lib @@ -0,0 +1,99 @@ +!!index array str + "CYS" +!entry.CYS.unit.atoms table str name str type int typex int resx int flags int seq int elmnt dbl chg + "C" "C" 0 1 131072 1 6 0.000000 + "O" "O" 0 1 131072 2 8 0.000000 + "N" "N" 0 1 131072 3 7 0.000000 + "H" "H" 0 1 131072 4 1 0.000000 + "CA" "C" 0 1 131072 5 6 0.000000 + "HA" "H" 0 1 131072 6 1 0.000000 + "CB" "C" 0 1 131072 7 6 0.000000 + "HB2" "H" 0 1 131072 8 1 0.000000 + "HB3" "H" 0 1 131072 9 1 0.000000 + "SG" "S" 0 1 131072 10 16 0.000000 + "HG" "H" 0 1 131072 11 1 0.000000 +!entry.CYS.unit.atomspertinfo table str pname str ptype int ptypex int pelmnt dbl pchg + "C" "C" 0 -1 0.0 + "O" "O" 0 -1 0.0 + "N" "N" 0 -1 0.0 + "H" "H" 0 -1 0.0 + "CA" "C" 0 -1 0.0 + "HA" "H" 0 -1 0.0 + "CB" "C" 0 -1 0.0 + "HB2" "H" 0 -1 0.0 + "HB3" "H" 0 -1 0.0 + "SG" "S" 0 -1 0.0 + "HG" "H" 0 -1 0.0 +!entry.CYS.unit.boundbox array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.CYS.unit.childsequence single int + 2 +!entry.CYS.unit.connect array int + 0 + 0 +!entry.CYS.unit.connectivity table int atom1x int atom2x int flags + 1 2 1 + 3 4 1 + 3 5 1 + 5 6 1 + 5 7 1 + 5 1 1 + 7 8 1 + 7 9 1 + 7 10 1 + 10 11 1 +!entry.CYS.unit.hierarchy table str abovetype int abovex str belowtype int belowx + "U" 0 "R" 1 + "R" 1 "A" 1 + "R" 1 "A" 2 + "R" 1 "A" 3 + "R" 1 "A" 4 + "R" 1 "A" 5 + "R" 1 "A" 6 + "R" 1 "A" 7 + "R" 1 "A" 8 + "R" 1 "A" 9 + "R" 1 "A" 10 + "R" 1 "A" 11 +!entry.CYS.unit.name single str + "CYS" +!entry.CYS.unit.positions table dbl x dbl y dbl z + 2.200000 1.200000 0.200000 + 3.400000 1.200000 0.400000 + 0.000000 0.000000 0.000000 + 0.000000 0.000000 1.010000 + 1.450000 0.000000 -0.500000 + 1.450000 0.890000 -1.090000 + 2.200000 -1.200000 0.200000 + 1.700000 -2.100000 -0.100000 + 3.200000 -1.100000 -0.200000 + 2.400000 -1.400000 2.000000 + 3.500000 -1.100000 2.300000 +!entry.CYS.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x + 0 0 0 0 0 0 +!entry.CYS.unit.residues table str name int seq int childseq int startatomx str restype int imagingx + "CYS" 1 12 1 "?" 0 +!entry.CYS.unit.residuesPdbSequenceNumber array int + 0 +!entry.CYS.unit.solventcap array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.CYS.unit.velocities table dbl x dbl y dbl z + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 diff --git a/test/Test_TIMap/daa_identity.map.save b/test/Test_TIMap/daa_identity.map.save new file mode 100644 index 0000000000..e1e6fce1d0 --- /dev/null +++ b/test/Test_TIMap/daa_identity.map.save @@ -0,0 +1,26 @@ +# timap tgt='DAA' template='DAA' +# kind tgt=base template=base +# mapped= 22 insertion= 0 unmapped_template= 0 n_tgt= 22 n_tpl= 22 +#Out TgtAt TgtName TplAt TplName Ins + 1 1 O4' 1 O4' 0 + 2 2 C1' 2 C1' 0 + 3 3 H1' 3 H1' 0 + 4 4 N9 4 N9 0 + 5 5 C8 5 C8 0 + 6 6 H8 6 H8 0 + 7 7 N7 7 N7 0 + 8 8 C5 8 C5 0 + 9 9 C6 9 C6 0 + 10 10 N6 10 N6 0 + 11 11 H61 11 H61 0 + 12 12 H62 12 H62 0 + 13 13 N1 13 N1 0 + 14 14 C2 14 C2 0 + 15 15 H2 15 H2 0 + 16 16 N3 16 N3 0 + 17 17 C4 17 C4 0 + 18 18 C2' 18 C2' 0 + 19 19 H2'1 19 H2'1 0 + 20 20 HHO4 20 HHO4 0 + 21 21 H2'2 21 H2'2 0 + 22 22 H2'3 22 H2'3 0 diff --git a/test/Test_TIMap/ern_naorder.map.save b/test/Test_TIMap/ern_naorder.map.save new file mode 100644 index 0000000000..90879ea665 --- /dev/null +++ b/test/Test_TIMap/ern_naorder.map.save @@ -0,0 +1,36 @@ +# timap tgt='ERN[ERN]' template='ERN[ERN]' +# kind tgt=nucleotide template=nucleotide +# mapped= 32 insertion= 0 unmapped_template= 0 n_tgt= 32 n_tpl= 32 +#Out TgtAt TgtName TplAt TplName Ins + 1 15 P 15 P 0 + 2 16 OP1 16 OP1 0 + 3 17 OP2 17 OP2 0 + 4 18 O5' 18 O5' 0 + 5 19 C5' 19 C5' 0 + 6 20 H5' 20 H5' 0 + 7 21 H5'' 21 H5'' 0 + 8 22 C4' 22 C4' 0 + 9 23 H4' 23 H4' 0 + 10 24 O4' 24 O4' 0 + 11 25 C1' 25 C1' 0 + 12 31 H1' 31 H1' 0 + 13 1 N9 1 N9 0 + 14 14 C4 14 C4 0 + 15 2 C8 2 C8 0 + 16 3 H8 3 H8 0 + 17 5 C5 5 C5 0 + 18 13 N3 13 N3 0 + 19 4 N7 4 N7 0 + 20 6 C6 6 C6 0 + 21 11 C2 11 C2 0 + 22 12 H2 12 H2 0 + 23 10 N1 10 N1 0 + 24 7 N6 7 N6 0 + 25 8 H61 8 H61 0 + 26 9 H62 9 H62 0 + 27 26 C3' 26 C3' 0 + 28 27 H3' 27 H3' 0 + 29 28 C2' 28 C2' 0 + 30 29 H2' 29 H2' 0 + 31 32 H2'' 32 H2'' 0 + 32 30 O3' 30 O3' 0 diff --git a/test/Test_TIMap/fle_ern.atoms.save b/test/Test_TIMap/fle_ern.atoms.save new file mode 100644 index 0000000000..b5726a215a --- /dev/null +++ b/test/Test_TIMap/fle_ern.atoms.save @@ -0,0 +1,37 @@ +# timap dual-topology atoms +# Kind: SHARED = real in both; TPL_ONLY = dummy in lambda 1; TGT_ONLY = dummy in lambda 0 +#Idx Kind Name0 Q0 Name1 Q1 + 1 SHARED N9 -0.026798 N9 -0.025018 + 2 SHARED C8 0.160711 C8 0.201257 + 3 SHARED H8 0.187713 H8 0.155808 + 4 SHARED N7 -0.617458 N7 -0.605312 + 5 SHARED C5 0.072505 C5 0.051669 + 6 SHARED C6 0.689747 C6 0.703195 + 7 SHARED N6 -0.912238 N6 -0.898947 + 8 SHARED H61 0.416728 H61 0.412847 + 9 SHARED H62 0.416728 H62 0.412847 + 10 SHARED N1 -0.762348 N1 -0.759007 + 11 SHARED C2 0.571639 C2 0.589424 + 12 SHARED H2 0.059804 H2 0.047455 + 13 SHARED N3 -0.741649 N3 -0.697409 + 14 SHARED C4 0.380026 C4 0.306300 + 15 SHARED P 1.166744 P 1.166744 + 16 SHARED OP1 -0.775638 OP1 -0.775638 + 17 SHARED OP2 -0.775638 OP2 -0.775638 + 18 SHARED O5' -0.498667 O5' -0.498667 + 19 SHARED C5' -0.006896 C5' 0.055388 + 20 SHARED H5' 0.075442 H5' 0.067399 + 21 SHARED H5'' 0.075442 H5'' 0.067399 + 22 SHARED C4' 0.162991 C4' 0.105714 + 23 SHARED H4' 0.117666 H4' 0.116533 + 24 SHARED O4' -0.368894 O4' -0.357420 + 25 SHARED C1' 0.043124 C1' 0.039109 + 26 SHARED C3' 0.071340 C3' 0.200707 + 27 SHARED H3' 0.098555 H3' 0.061046 + 28 SHARED C2' -0.085352 C2' 0.066505 + 29 TGT_ONLY O2' 0.000000 O2' -0.618434 + 30 TGT_ONLY HO2' 0.000000 HO2' 0.415509 + 31 SHARED H2' 0.071840 H2' 0.096482 + 32 SHARED O3' -0.522908 O3' -0.527064 + 33 SHARED H1' 0.183902 H1' 0.199218 + 34 TPL_ONLY H2'' 0.071840 H2'' 0.000000 diff --git a/test/Test_TIMap/fle_ern.scmask.save b/test/Test_TIMap/fle_ern.scmask.save new file mode 100644 index 0000000000..fb88b72c41 --- /dev/null +++ b/test/Test_TIMap/fle_ern.scmask.save @@ -0,0 +1,7 @@ +# timap dual-topology TI masks +# n_dual= 34 n_shared= 31 dummy_in_lambda0= 2 dummy_in_lambda1= 1 +# lambda 0 = template (real unmatched, dummy insertions) +# lambda 1 = target (real insertions, dummy unmatched) +# pmemd: scmask1 is unique to lambda 0; scmask2 is unique to lambda 1. +scmask1=":1@H2''" +scmask2=":1@O2',:1@HO2'" diff --git a/test/Test_TIMap/fle_to_ern.map.save b/test/Test_TIMap/fle_to_ern.map.save new file mode 100644 index 0000000000..0589b9d6f6 --- /dev/null +++ b/test/Test_TIMap/fle_to_ern.map.save @@ -0,0 +1,39 @@ +# timap tgt='FLE[FLE]' template='ERN[ERN]' +# kind tgt=nucleotide template=nucleotide +# mapped= 31 insertion= 2 unmapped_template= 1 n_tgt= 33 n_tpl= 32 +#Out TgtAt TgtName TplAt TplName Ins + 1 1 N9 1 N9 0 + 2 2 C8 2 C8 0 + 3 3 H8 3 H8 0 + 4 4 N7 4 N7 0 + 5 5 C5 5 C5 0 + 6 6 C6 6 C6 0 + 7 7 N6 7 N6 0 + 8 8 H61 8 H61 0 + 9 9 H62 9 H62 0 + 10 10 N1 10 N1 0 + 11 11 C2 11 C2 0 + 12 12 H2 12 H2 0 + 13 13 N3 13 N3 0 + 14 14 C4 14 C4 0 + 15 15 P 15 P 0 + 16 16 OP1 16 OP1 0 + 17 17 OP2 17 OP2 0 + 18 18 O5' 18 O5' 0 + 19 19 C5' 19 C5' 0 + 20 20 H5' 20 H5' 0 + 21 21 H5'' 21 H5'' 0 + 22 22 C4' 22 C4' 0 + 23 23 H4' 23 H4' 0 + 24 24 O4' 24 O4' 0 + 25 25 C1' 25 C1' 0 + 26 26 C3' 26 C3' 0 + 27 27 H3' 27 H3' 0 + 28 28 C2' 28 C2' 0 + 29 30 O2' --- --- 1 + 30 33 HO2' --- --- 1 + 31 29 H2' 29 H2' 0 + 32 31 O3' 30 O3' 0 + 33 32 H1' 31 H1' 0 +# unmapped template atoms: +# 32 H2'' diff --git a/test/Test_TIMap/path_mix.map.save b/test/Test_TIMap/path_mix.map.save new file mode 100644 index 0000000000..8ca363362a --- /dev/null +++ b/test/Test_TIMap/path_mix.map.save @@ -0,0 +1,19 @@ +# timap tgt='phenol.mol2' template='benzene.pdb' +# kind tgt=unknown template=unknown +# mapped= 11 insertion= 2 unmapped_template= 1 n_tgt= 13 n_tpl= 12 +#Out TgtAt TgtName TplAt TplName Ins + 1 1 C1 1 C1 0 + 2 2 O1 --- --- 1 + 3 3 HO1 --- --- 1 + 4 4 C2 3 C2 0 + 5 5 H2 4 H2 0 + 6 6 C3 5 C3 0 + 7 7 H3 6 H3 0 + 8 8 C4 7 C4 0 + 9 9 H4 8 H4 0 + 10 10 C5 9 C5 0 + 11 11 H5 10 H5 0 + 12 12 C6 11 C6 0 + 13 13 H6 12 H6 0 +# unmapped template atoms: +# 2 H1 diff --git a/test/Test_TIMap/phenol.mol2 b/test/Test_TIMap/phenol.mol2 new file mode 100644 index 0000000000..1e3eb200ee --- /dev/null +++ b/test/Test_TIMap/phenol.mol2 @@ -0,0 +1,37 @@ +@MOLECULE +phenol analog of benzene (OH replaces H1) + 13 13 1 0 0 +SMALL +NO_CHARGES + + +@ATOM + 1 C1 1.2000 0.0000 0.0000 C 1 PHE 0.000000 + 2 O1 2.4000 0.0000 0.0000 O 1 PHE 0.000000 + 3 HO1 2.7600 0.8900 0.0000 H 1 PHE 0.000000 + 4 C2 0.6000 1.0390 0.0000 C 1 PHE 0.000000 + 5 H2 1.1400 1.9740 0.0000 H 1 PHE 0.000000 + 6 C3 -0.6000 1.0390 0.0000 C 1 PHE 0.000000 + 7 H3 -1.1400 1.9740 0.0000 H 1 PHE 0.000000 + 8 C4 -1.2000 0.0000 0.0000 C 1 PHE 0.000000 + 9 H4 -2.2800 0.0000 0.0000 H 1 PHE 0.000000 + 10 C5 -0.6000 -1.0390 0.0000 C 1 PHE 0.000000 + 11 H5 -1.1400 -1.9740 0.0000 H 1 PHE 0.000000 + 12 C6 0.6000 -1.0390 0.0000 C 1 PHE 0.000000 + 13 H6 1.1400 -1.9740 0.0000 H 1 PHE 0.000000 +@BOND + 1 1 2 1 + 2 2 3 1 + 3 1 4 1 + 4 4 5 1 + 5 4 6 1 + 6 6 7 1 + 7 6 8 1 + 8 8 9 1 + 9 8 10 1 + 10 10 11 1 + 11 10 12 1 + 12 12 13 1 + 13 12 1 1 +@SUBSTRUCTURE + 1 PHE 1 **** 0 **** **** diff --git a/test/Test_TIMap/phenol.pdb b/test/Test_TIMap/phenol.pdb new file mode 100644 index 0000000000..a115699478 --- /dev/null +++ b/test/Test_TIMap/phenol.pdb @@ -0,0 +1,27 @@ +ATOM 1 C1 PHE A 1 1.200 0.000 0.000 1.00 0.00 C +ATOM 2 O1 PHE A 1 2.400 0.000 0.000 1.00 0.00 O +ATOM 3 HO1 PHE A 1 2.760 0.890 0.000 1.00 0.00 H +ATOM 4 C2 PHE A 1 0.600 1.039 0.000 1.00 0.00 C +ATOM 5 H2 PHE A 1 1.140 1.974 0.000 1.00 0.00 H +ATOM 6 C3 PHE A 1 -0.600 1.039 0.000 1.00 0.00 C +ATOM 7 H3 PHE A 1 -1.140 1.974 0.000 1.00 0.00 H +ATOM 8 C4 PHE A 1 -1.200 0.000 0.000 1.00 0.00 C +ATOM 9 H4 PHE A 1 -2.280 0.000 0.000 1.00 0.00 H +ATOM 10 C5 PHE A 1 -0.600 -1.039 0.000 1.00 0.00 C +ATOM 11 H5 PHE A 1 -1.140 -1.974 0.000 1.00 0.00 H +ATOM 12 C6 PHE A 1 0.600 -1.039 0.000 1.00 0.00 C +ATOM 13 H6 PHE A 1 1.140 -1.974 0.000 1.00 0.00 H +CONECT 1 2 +CONECT 2 3 +CONECT 1 4 +CONECT 4 5 +CONECT 4 6 +CONECT 6 7 +CONECT 6 8 +CONECT 8 9 +CONECT 8 10 +CONECT 10 11 +CONECT 10 12 +CONECT 12 13 +CONECT 12 1 +END diff --git a/test/Test_TIMap/phenol_bnz.atoms.save b/test/Test_TIMap/phenol_bnz.atoms.save new file mode 100644 index 0000000000..15c028c0b6 --- /dev/null +++ b/test/Test_TIMap/phenol_bnz.atoms.save @@ -0,0 +1,17 @@ +# timap dual-topology atoms +# Kind: SHARED = real in both; TPL_ONLY = dummy in lambda 1; TGT_ONLY = dummy in lambda 0 +#Idx Kind Name0 Q0 Name1 Q1 + 1 SHARED C1 0.000000 C1 0.000000 + 2 TGT_ONLY O1 0.000000 O1 0.000000 + 3 TGT_ONLY HO1 0.000000 HO1 0.000000 + 4 TPL_ONLY H1 0.000000 H1 0.000000 + 5 SHARED C2 0.000000 C2 0.000000 + 6 SHARED H2 0.000000 H2 0.000000 + 7 SHARED C3 0.000000 C3 0.000000 + 8 SHARED H3 0.000000 H3 0.000000 + 9 SHARED C4 0.000000 C4 0.000000 + 10 SHARED H4 0.000000 H4 0.000000 + 11 SHARED C5 0.000000 C5 0.000000 + 12 SHARED H5 0.000000 H5 0.000000 + 13 SHARED C6 0.000000 C6 0.000000 + 14 SHARED H6 0.000000 H6 0.000000 diff --git a/test/Test_TIMap/phenol_bnz.scmask.save b/test/Test_TIMap/phenol_bnz.scmask.save new file mode 100644 index 0000000000..ceb1365ef9 --- /dev/null +++ b/test/Test_TIMap/phenol_bnz.scmask.save @@ -0,0 +1,7 @@ +# timap dual-topology TI masks +# n_dual= 14 n_shared= 11 dummy_in_lambda0= 2 dummy_in_lambda1= 1 +# lambda 0 = template (real unmatched, dummy insertions) +# lambda 1 = target (real insertions, dummy unmatched) +# pmemd: scmask1 is unique to lambda 0; scmask2 is unique to lambda 1. +scmask1=":1@H1" +scmask2=":1@O1,:1@HO1" diff --git a/test/Test_TIMap/phenol_to_benzene.map.save b/test/Test_TIMap/phenol_to_benzene.map.save new file mode 100644 index 0000000000..c43c925c56 --- /dev/null +++ b/test/Test_TIMap/phenol_to_benzene.map.save @@ -0,0 +1,19 @@ +# timap tgt='phenol' template='benzene' +# kind tgt=unknown template=unknown +# mapped= 11 insertion= 2 unmapped_template= 1 n_tgt= 13 n_tpl= 12 +#Out TgtAt TgtName TplAt TplName Ins + 1 1 C1 1 C1 0 + 2 2 O1 --- --- 1 + 3 3 HO1 --- --- 1 + 4 4 C2 3 C2 0 + 5 5 H2 4 H2 0 + 6 6 C3 5 C3 0 + 7 7 H3 6 H3 0 + 8 8 C4 7 C4 0 + 9 9 H4 8 H4 0 + 10 10 C5 9 C5 0 + 11 11 H5 10 H5 0 + 12 12 C6 11 C6 0 + 13 13 H6 12 H6 0 +# unmapped template atoms: +# 2 H1 diff --git a/test/Test_TIMap/sec.lib b/test/Test_TIMap/sec.lib new file mode 100644 index 0000000000..951f53dc09 --- /dev/null +++ b/test/Test_TIMap/sec.lib @@ -0,0 +1,99 @@ +!!index array str + "SEC" +!entry.SEC.unit.atoms table str name str type int typex int resx int flags int seq int elmnt dbl chg + "N" "N" 0 1 131072 1 7 0.000000 + "H" "H" 0 1 131072 2 1 0.000000 + "CA" "C" 0 1 131072 3 6 0.000000 + "HA" "H" 0 1 131072 4 1 0.000000 + "CB" "C" 0 1 131072 5 6 0.000000 + "HB2" "H" 0 1 131072 6 1 0.000000 + "HB3" "H" 0 1 131072 7 1 0.000000 + "SE" "Se" 0 1 131072 8 34 0.000000 + "HG" "H" 0 1 131072 9 1 0.000000 + "C" "C" 0 1 131072 10 6 0.000000 + "O" "O" 0 1 131072 11 8 0.000000 +!entry.SEC.unit.atomspertinfo table str pname str ptype int ptypex int pelmnt dbl pchg + "N" "N" 0 -1 0.0 + "H" "H" 0 -1 0.0 + "CA" "C" 0 -1 0.0 + "HA" "H" 0 -1 0.0 + "CB" "C" 0 -1 0.0 + "HB2" "H" 0 -1 0.0 + "HB3" "H" 0 -1 0.0 + "SE" "Se" 0 -1 0.0 + "HG" "H" 0 -1 0.0 + "C" "C" 0 -1 0.0 + "O" "O" 0 -1 0.0 +!entry.SEC.unit.boundbox array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.SEC.unit.childsequence single int + 2 +!entry.SEC.unit.connect array int + 0 + 0 +!entry.SEC.unit.connectivity table int atom1x int atom2x int flags + 1 2 1 + 1 3 1 + 3 4 1 + 3 5 1 + 3 10 1 + 5 6 1 + 5 7 1 + 5 8 1 + 8 9 1 + 10 11 1 +!entry.SEC.unit.hierarchy table str abovetype int abovex str belowtype int belowx + "U" 0 "R" 1 + "R" 1 "A" 1 + "R" 1 "A" 2 + "R" 1 "A" 3 + "R" 1 "A" 4 + "R" 1 "A" 5 + "R" 1 "A" 6 + "R" 1 "A" 7 + "R" 1 "A" 8 + "R" 1 "A" 9 + "R" 1 "A" 10 + "R" 1 "A" 11 +!entry.SEC.unit.name single str + "SEC" +!entry.SEC.unit.positions table dbl x dbl y dbl z + 0.000000 0.000000 0.000000 + 0.000000 0.000000 1.010000 + 1.450000 0.000000 -0.500000 + 1.450000 0.890000 -1.090000 + 2.200000 -1.200000 0.200000 + 1.700000 -2.100000 -0.100000 + 3.200000 -1.100000 -0.200000 + 2.500000 -1.500000 2.200000 + 3.700000 -1.100000 2.500000 + 2.200000 1.200000 0.200000 + 3.400000 1.200000 0.400000 +!entry.SEC.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x + 0 0 0 0 0 0 +!entry.SEC.unit.residues table str name int seq int childseq int startatomx str restype int imagingx + "SEC" 1 12 1 "?" 0 +!entry.SEC.unit.residuesPdbSequenceNumber array int + 0 +!entry.SEC.unit.solventcap array dbl + -1.000000 + 0.0 + 0.0 + 0.0 + 0.0 +!entry.SEC.unit.velocities table dbl x dbl y dbl z + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 diff --git a/test/Test_TIMap/sec.mol2 b/test/Test_TIMap/sec.mol2 new file mode 100644 index 0000000000..07cc7f92f7 --- /dev/null +++ b/test/Test_TIMap/sec.mol2 @@ -0,0 +1,32 @@ +@MOLECULE +selenocysteine analog of cysteine + 11 10 1 0 0 +SMALL +NO_CHARGES + + +@ATOM + 1 N 0.0000 0.0000 0.0000 N 1 SEC 0.000000 + 2 H 0.0000 0.0000 1.0100 H 1 SEC 0.000000 + 3 CA 1.4500 0.0000 -0.5000 C 1 SEC 0.000000 + 4 HA 1.4500 0.8900 -1.0900 H 1 SEC 0.000000 + 5 CB 2.2000 -1.2000 0.2000 C 1 SEC 0.000000 + 6 HB2 1.7000 -2.1000 -0.1000 H 1 SEC 0.000000 + 7 HB3 3.2000 -1.1000 -0.2000 H 1 SEC 0.000000 + 8 SE 2.5000 -1.5000 2.2000 Se 1 SEC 0.000000 + 9 HG 3.7000 -1.1000 2.5000 H 1 SEC 0.000000 + 10 C 2.2000 1.2000 0.2000 C 1 SEC 0.000000 + 11 O 3.4000 1.2000 0.4000 O 1 SEC 0.000000 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 3 4 1 + 4 3 5 1 + 5 3 10 1 + 6 5 6 1 + 7 5 7 1 + 8 5 8 1 + 9 8 9 1 + 10 10 11 1 +@SUBSTRUCTURE + 1 SEC 1 **** 0 **** **** diff --git a/test/Test_TIMap/sec_cys.atoms.save b/test/Test_TIMap/sec_cys.atoms.save new file mode 100644 index 0000000000..7983b5ec5e --- /dev/null +++ b/test/Test_TIMap/sec_cys.atoms.save @@ -0,0 +1,15 @@ +# timap dual-topology atoms +# Kind: SHARED = real in both; TPL_ONLY = dummy in lambda 1; TGT_ONLY = dummy in lambda 0 +#Idx Kind Name0 Q0 Name1 Q1 + 1 SHARED N 0.000000 N 0.000000 + 2 SHARED H 0.000000 H 0.000000 + 3 SHARED CA 0.000000 CA 0.000000 + 4 SHARED HA 0.000000 HA 0.000000 + 5 SHARED CB 0.000000 CB 0.000000 + 6 TGT_ONLY SE 0.000000 SE 0.000000 + 7 SHARED HB2 0.000000 HB2 0.000000 + 8 SHARED HB3 0.000000 HB3 0.000000 + 9 TPL_ONLY SG 0.000000 SG 0.000000 + 10 SHARED HG 0.000000 HG 0.000000 + 11 SHARED C 0.000000 C 0.000000 + 12 SHARED O 0.000000 O 0.000000 diff --git a/test/Test_TIMap/sec_cys.scmask.save b/test/Test_TIMap/sec_cys.scmask.save new file mode 100644 index 0000000000..d93870f150 --- /dev/null +++ b/test/Test_TIMap/sec_cys.scmask.save @@ -0,0 +1,7 @@ +# timap dual-topology TI masks +# n_dual= 12 n_shared= 10 dummy_in_lambda0= 1 dummy_in_lambda1= 1 +# lambda 0 = template (real unmatched, dummy insertions) +# lambda 1 = target (real insertions, dummy unmatched) +# pmemd: scmask1 is unique to lambda 0; scmask2 is unique to lambda 1. +scmask1=":1@SG" +scmask2=":1@SE" diff --git a/test/Test_TIMap/sec_to_cys.map.save b/test/Test_TIMap/sec_to_cys.map.save new file mode 100644 index 0000000000..a2d8ac2bcf --- /dev/null +++ b/test/Test_TIMap/sec_to_cys.map.save @@ -0,0 +1,17 @@ +# timap tgt='SEC[SEC]' template='CYS[CYS]' +# kind tgt=amino template=amino +# mapped= 10 insertion= 1 unmapped_template= 1 n_tgt= 11 n_tpl= 11 +#Out TgtAt TgtName TplAt TplName Ins + 1 1 N 1 N 0 + 2 2 H 2 H 0 + 3 3 CA 3 CA 0 + 4 4 HA 4 HA 0 + 5 5 CB 5 CB 0 + 6 8 SE --- --- 1 + 7 6 HB2 6 HB2 0 + 8 7 HB3 7 HB3 0 + 9 9 HG 9 HG 0 + 10 10 C 10 C 0 + 11 11 O 11 O 0 +# unmapped template atoms: +# 8 SG diff --git a/test/Test_TIMap/series_map_phenol.map.save b/test/Test_TIMap/series_map_phenol.map.save new file mode 100644 index 0000000000..c43c925c56 --- /dev/null +++ b/test/Test_TIMap/series_map_phenol.map.save @@ -0,0 +1,19 @@ +# timap tgt='phenol' template='benzene' +# kind tgt=unknown template=unknown +# mapped= 11 insertion= 2 unmapped_template= 1 n_tgt= 13 n_tpl= 12 +#Out TgtAt TgtName TplAt TplName Ins + 1 1 C1 1 C1 0 + 2 2 O1 --- --- 1 + 3 3 HO1 --- --- 1 + 4 4 C2 3 C2 0 + 5 5 H2 4 H2 0 + 6 6 C3 5 C3 0 + 7 7 H3 6 H3 0 + 8 8 C4 7 C4 0 + 9 9 H4 8 H4 0 + 10 10 C5 9 C5 0 + 11 11 H5 10 H5 0 + 12 12 C6 11 C6 0 + 13 13 H6 12 H6 0 +# unmapped template atoms: +# 2 H1 diff --git a/test/Test_TIMap/series_parent.dat.save b/test/Test_TIMap/series_parent.dat.save new file mode 100644 index 0000000000..f0ab9afd7b --- /dev/null +++ b/test/Test_TIMap/series_parent.dat.save @@ -0,0 +1,7 @@ +# timap series parent selection +# score = sum_j exp(-0.1*(insertions+unmapped)) mapping j onto candidate +# adapted from util/leadopt/similarity.py (Redmann & Summa / Summa Lab) +parent benzene +#Idx Name Natom SimSum MappedSum + 1 phenol 13 0.740818 11 + 2 benzene 12 0.740818 11 parent