/* FILE: base_information.C */ /* VERSION: 1.600 */ /* DATE: 2/3/2005 */ /* Comments and Questions to: sw-help@rcsb.rutgers.edu */ /* COPYRIGHT 1999-2005 Rutgers - The State University of New Jersey This software is provided WITHOUT WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE OR ANY OTHER WARRANTY, EXPRESS OR IMPLIED. RUTGERS MAKE NO REPRESENTATION OR WARRANTY THAT THE SOFTWARE WILL NOT INFRINGE ANY PATENT, COPYRIGHT OR OTHER PROPRIETARY RIGHT. The user of this software shall indemnify, hold harmless and defend Rutgers, its governors, trustees, officers, employees, students, agents and the authors against any and all claims, suits, losses, liabilities, damages, costs, fees, and expenses including reasonable attorneys' fees resulting from or arising out of the use of this software. This indemnification shall include, but is not limited to, any and all claims alleging products liability. */ /* PDB SOFTWARE LICENSE AGREEMENT BY CLICKING THE ACCEPTANCE BUTTON OR INSTALLING OR USING THIS "SOFTWARE, THE INDIVIDUAL OR ENTITY LICENSING THE SOFTWARE ("LICENSEE") IS CONSENTING TO BE BOUND BY AND IS BECOMING A PARTY TO THIS AGREEMENT. IF LICENSEE DOES NOT AGREE TO ALL OF THE TERMS OF THIS AGREEMENT THE LICENSEE MUST NOT INSTALL OR USE THE SOFTWARE. 1. LICENSE AGREEMENT This is a license between you ("Licensee") and the Protein Data Bank (PDB) at Rutgers, The State University of New Jersey (hereafter referred to as "RUTGERS"). The software is owned by RUTGERS and protected by copyright laws, and some elements are protected by laws governing trademarks, trade dress and trade secrets, and may be protected by patent laws. 2. LICENSE GRANT RUTGERS grants you, and you hereby accept, non-exclusive, royalty-free perpetual license to install, use, modify, prepare derivative works, incorporate into other computer software, and distribute in binary and source code format, or any derivative work thereof, together with any associated media, printed materials, and on-line or electronic documentation (if any) provided by RUTGERS (collectively, the "SOFTWARE"), subject to the following terms and conditions: (i) any distribution of the SOFTWARE shall bind the receiver to the terms and conditions of this Agreement; (ii) any distribution of the SOFTWARE in modified form shall clearly state that the SOFTWARE has been modified from the version originally obtained from RUTGERS. 2. COPYRIGHT; RETENTION OF RIGHTS. The above license grant is conditioned on the following: (i) you must reproduce all copyright notices and other proprietary notices on any copies of the SOFTWARE and you must not remove such notices; (ii) in the event you compile the SOFTWARE, you will include the copyright notice with the binary in such a manner as to allow it to be easily viewable; (iii) if you incorporate the SOFTWARE into other code, you must provide notice that the code contains the SOFTWARE and include a copy of the copyright notices and other proprietary notices. All copies of the SOFTWARE shall be subject to the terms of this Agreement. 3. NO MAINTENANCE OR SUPPORT; TREATMENT OF ENHANCEMENTS RUTGERS is under no obligation whatsoever to: (i) provide maintenance or support for the SOFTWARE; or (ii) to notify you of bug fixes, patches, or upgrades to the features, functionality or performance of the SOFTWARE ("Enhancements") (if any), whether developed by RUTGERS or third parties. If, in its sole discretion, RUTGERS makes an Enhancement available to you and RUTGERS does not separately enter into a written license agreement with you relating to such bug fix, patch or upgrade, then it shall be deemed incorporated into the SOFTWARE and subject to this Agreement. You are under no obligation whatsoever to provide any Enhancements to RUTGERS or the public that you may develop over time; however, if you choose to provide your Enhancements to RUTGERS, or if you choose to otherwise publish or distribute your Enhancements, in source code form without contemporaneously requiring end users or RUTGERS to enter into a separate written license agreement for such Enhancements, then you hereby grant RUTGERS a non-exclusive, royalty-free perpetual license to install, use, modify, prepare derivative works, incorporate into the SOFTWARE or other computer software, distribute, and sublicense your Enhancements or derivative works thereof, in binary and source code form. 4. FEES. There is no license fee for the SOFTWARE. If Licensee wishes to receive the SOFTWARE on media, there may be a small charge for the media and for shipping and handling. Licensee is responsible for any and all taxes. 5. TERMINATION. Without prejudice to any other rights, Licensor may terminate this Agreement if Licensee breaches any of its terms and conditions. Upon termination, Licensee shall destroy all copies of the SOFTWARE. 6. PROPRIETARY RIGHTS. Title, ownership rights, and intellectual property rights in the Product shall remain with RUTGERS. Licensee acknowledges such ownership and intellectual property rights and will not take any action to jeopardize, limit or interfere in any manner with RUTGERS' ownership of or rights with respect to the SOFTWARE. The SOFTWARE is protected by copyright and other intellectual property laws and by international treaties. Title and related rights in the content accessed through the SOFTWARE is the property of the applicable content owner and is protected by applicable law. The license granted under this Agreement gives Licensee no rights to such content. 7. DISCLAIMER OF WARRANTY. THE SOFTWARE IS PROVIDED FREE OF CHARGE, AND, THEREFORE, ON AN "AS IS" BASIS, WITHOUT WARRANTY OF ANY KIND, INCLUDING WITHOUT LIMITATION THE WARRANTIES THAT IT IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE SOFTWARE IS BORNE BY LICENSEE. SHOULD THE SOFTWARE PROVE DEFECTIVE IN ANY RESPECT, THE LICENSEE AND NOT LICENSOR ASSUMES THE ENTIRE COST OF ANY SERVICE AND REPAIR. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL PART OF THIS AGREEMENT. NO USE OF THE PRODUCT IS AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER. 8. LIMITATION OF LIABILITY. TO THE MAXIMUM EXTENT PERMITTED BY APPLICABLE LAW, IN NO EVENT WILL LICENSOR BE LIABLE FOR ANY INDIRECT, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF OR INABILITY TO USE THE SOFTWARE, INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK STOPPAGE, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF ADVISED OF THE POSSIBILITY THEREOF. */ #include #include #include #include #include #include #include "_bpindex.h" #include "_index.h" #include "BTreeObj.h" #include "TblIndexObj.h" #include "ndb_misclib.h" #include "_validation.h" #include "base_information.h" // #include "contact.h" #include "misc_util.h" #include "xtal.h" #include "util.h" static void realloc_base(_BASE*** _base, const int n_base, int &alloc_base); static _BASE *base_frame(_RESIDUE *res, const int op, const int lx, const int ly, const int lz, const NDBCELLCONST *cell); static int find_match(_RESIDUE *residue, double *eRing_xyz, double *sRing_xyz, double **block, int& Y_or_R); static _BASE_PAIR *find_all_pairs(const int n_base, _BASE **base); static CONTACT **find_real_base_pairs(const int n_base, _BASE **base, _BASE_PAIR *base_pair, int &n_contact, const NDBCELLCONST *cell, _MOLECULE *mol); static double cal_dist(_ATOM *atom_a, _ATOM *atom_b, const NDBCELLCONST *cell, const int op, const int lx, const int ly, const int lz); static int get_best_pair(_BASE_PAIR *base_pair, int &pair_type, int &type_i, int &type_j, int &cis_or_trans, int *skip); static void bpstep_par(COORD rot1[3], const COORD org1, COORD rot2[3], const COORD org2, double *pars, COORD mst_orien[3], COORD &mst_org); static void helical_par(COORD rot1[3], const COORD org1, COORD rot2[3], const COORD org2, double *pars, COORD mst_orien[3], COORD &mst_org); static void project_xyzP(const COORD P_1_i_plus_1, const COORD P_2_i, COORD mst_orien[3], COORD mst_org, COORD& aveP); static double distance_ab(_BASE *a, const int ia, _BASE *b, const int ib); static void get_ij(const int m, const int index, int *best_pairs, int &i, int &j); static void reverse(const int st, const int n, int *vec); static int is_linked(_BASE *base1, _BASE *base2); static int is_wc_geometry(_BASE *base1, _BASE *base2); static int check_others(_BASE **_base, const int i1, const int j1, const int i2, const int j2); static int check_direction(int *helix_idx, int *bp_idx, int *index, int *best_pairs, _BASE **_base, int *direction); static int find_atom_index(const int num, _ATOMNAME *name_list, const char *name); #ifdef SUN_OS #ifdef __cplusplus extern "C" { #endif #endif static int base_hbond_type_compare(const void *base_type1, const void *base_type2); static int name_compare1(const void *name1, const void *name2); static int base_type_compare(const void *base_type1, const void *base_type2); static int atom_edge_compare(const void *atom_edge1, const void *atom_edge2); #ifdef SUN_OS #ifdef __cplusplus } #endif #endif void _VALIDATION::get_bspair_information() { int i = 0, j = 0, max_count = 0, n_base = 0, alloc_base = 0, all_type[32]; int n_contact = 0, n_sym_op = 0, ret = -1, found, atomn_chain; xtal mycrys; _BASE **_base = NULL, *base = NULL; _CHAIN *chain = NULL; SYM_OP sym_op[100]; BTreeObj<_SIndex> _SymmetryIndex; _SIndex index; char buffer[100]; CONTACT **contact = NULL; if (_num_mol == 0) return; if (_cell == NULL) return; ; if (_bspair || _n_bspair) return; _classification = 0; atomn_chain = 0; chain = _molecules[0]->GetFirstChain(); while (chain) { if (chain->chain_type() == ATOMN_TOKEN) { atomn_chain++; _RESIDUE *res = chain->GetFirstResidue(); while (res) { if (!strcmp(res->alt_loc(), "") || res->alt_loc()[0] == 'A' || res->alt_loc()[0] == 'a' || res->alt_loc()[0] == '1') { base = base_frame(res, 0, 0, 0, 0, _cell); if (base) { realloc_base(&_base, n_base, alloc_base); _base[n_base] = base; n_base++; } } res = chain->GetNextResidue(); } } chain = _molecules[0]->GetNextChain(); } if (!n_base) return; for (i = 0; i < 32; i++) all_type[i] = 0; found = get_bspair_information(n_base, _base, all_type); for (i = 0; i < n_base; i++) { if (_base[i]) delete _base[i]; } delete [] _base; _base = NULL; n_base = 0; alloc_base = 0; if (atomn_chain == 1 && found) return; if (strcmp(_cell->space_group, "P 1")) { mycrys.setorth('c'); if (!mycrys.init(_cell, _molecules[0], 1)) return; mycrys.getSymmAtoms(); mycrys.moveAllAtomsToCenter(); mycrys.placeAtomsInGrid(); n_contact = 0; contact = get_contact(mycrys, n_contact , 2.2, 3.5, "sym", "kind atomn", "kind atomn", 0); if (contact) { n_sym_op = 0; ret = -1; for (i = 0; i < 100; i++) { sym_op[i].chn_id[0] = '\0'; sym_op[i].sym = 0; sym_op[i].transl[0] = 0; sym_op[i].transl[1] = 0; sym_op[i].transl[2] = 0; sym_op[i].count = 0; } for (i = 0; i < n_contact; i++) { if (contact[i]->sym <= 1 || strcmp(contact[i]->a_chnid, contact[i]->b_chnid)) continue; sprintf(buffer, "%s%d_%d_%d_%d", contact[i]->a_chnid, contact[i]->sym, contact[i]->transl[0], contact[i]->transl[1], contact[i]->transl[2]); index.set_value(buffer); ret = _SymmetryIndex.Seek(index); if (ret >= 0) { index = _SymmetryIndex.GetCurrData(); sym_op[index.pos()].count++; } else { index.set_pos(n_sym_op); _SymmetryIndex.Add(index); strcpy(sym_op[n_sym_op].chn_id, contact[i]->a_chnid); sym_op[n_sym_op].sym = contact[i]->sym - 1; sym_op[n_sym_op].transl[0] = contact[i]->transl[0]; sym_op[n_sym_op].transl[1] = contact[i]->transl[1]; sym_op[n_sym_op].transl[2] = contact[i]->transl[2]; sym_op[n_sym_op].count++; n_sym_op++; } } for (i = 0; i < n_contact; i++) if (contact[i]) delete contact[i]; delete [] contact; contact = NULL; j = 0; max_count = 0; for (i = 0; i < n_sym_op; i++) { if (sym_op[i].count > max_count) { max_count = sym_op[i].count; j = i; } } ret = 0; if (atomn_chain == 1 || (max_count > 20 && max_count > n_contact / 3)) { ret = 1; for (i = 0; i < _n_bb_params; i++) { // if (!strcmp(sym_op[j].chn_id, _bb_params[i]->ChnID_I) && // !strcmp(sym_op[j].chn_id, _bb_params[i]->ChnID_J)) { if (_bb_params[i]->ChnID_I == sym_op[j].chn_id && _bb_params[i]->ChnID_J == sym_op[j].chn_id) { ret = 0; break; } } } if (ret) { chain = _molecules[0]->GetFirstChain(); while (chain) { if (chain->chain_type() == ATOMN_TOKEN && !strcmp(chain->ChainID(), sym_op[j].chn_id)) { _RESIDUE *res = chain->GetFirstResidue(); while (res) { base = base_frame(res, 0, 0, 0, 0, _cell); if (base) { realloc_base(&_base, n_base, alloc_base); _base[n_base] = base; n_base++; } base = base_frame(res, sym_op[j].sym, sym_op[j].transl[0], sym_op[j].transl[1], sym_op[j].transl[2], _cell); if (base) { realloc_base(&_base, n_base, alloc_base); _base[n_base] = base; n_base++; } res = chain->GetNextResidue(); } } chain = _molecules[0]->GetNextChain(); } get_bspair_information(n_base, _base, all_type); for (i = 0; i < n_base; i++) { if (_base[i]) delete _base[i]; } delete [] _base; _base = NULL; n_base = 0; alloc_base = 0; } } } } int _VALIDATION::get_bspair_information(const int n_base, _BASE **_base, int *all_types) { int i = 0, j = 0, k = 0, l = 0, n_contact = 0, num_pair = 0, *index = NULL; int pair_type = -1, type_i, type_j, cis_or_trans, *best_pairs = NULL; _PAIR_TO **tmp = NULL; _BASE_PAIR *base_pair = find_all_pairs(n_base, _base); if (!base_pair) return 0; CONTACT **contact = find_real_base_pairs(n_base, _base, base_pair, n_contact, _cell, _molecules[0]); if (contact) { insert_base_pairs(n_contact, contact); for (i = 0; i < n_contact; i++) if (contact[i]) delete contact[i]; delete [] contact; contact = NULL; n_contact = 0; } for (i = 0; i < n_base; i++) { for (j = 0; j < base_pair[i].num; j++) { all_types[30]++; if (base_pair[i].pair_to[j]->type >= 0 && base_pair[i].pair_to[j]->type < 30) { all_types[31]++; all_types[base_pair[i].pair_to[j]->type]++; } } } for (i = 0; i < n_base - 1; i++) { for (k = 0; k < base_pair[i].num; k++) { if (base_pair[i].pair_to[k]->type < 0) continue; j = base_pair[i].pair_to[k]->index; for (l = 0; l < base_pair[j].num; l++) { if (base_pair[j].pair_to[l]->index == i) break; } if (l < base_pair[j].num) continue; if (base_pair[j].num >= base_pair[j].alloc_num) { base_pair[j].alloc_num += 2; tmp = base_pair[j].pair_to; base_pair[j].pair_to = new _PAIR_TO*[base_pair[j].alloc_num]; for (l = 0; l < base_pair[j].num; l++) base_pair[j].pair_to[l] = tmp[l]; for (l = base_pair[j].num; l < base_pair[j].alloc_num; l++) base_pair[j].pair_to[l] = NULL; delete tmp; } base_pair[j].pair_to[base_pair[j].num] = new _PAIR_TO; base_pair[j].pair_to[base_pair[j].num]->index = i; base_pair[j].pair_to[base_pair[j].num]->d2dp = base_pair[i].pair_to[k]->d2dp; base_pair[j].pair_to[base_pair[j].num]->type = base_pair[i].pair_to[k]->type; base_pair[j].pair_to[base_pair[j].num]->cis_or_trans = base_pair[i].pair_to[k]->cis_or_trans; base_pair[j].pair_to[base_pair[j].num]->type_i = base_pair[i].pair_to[k]->type_j; base_pair[j].pair_to[base_pair[j].num]->type_j = base_pair[i].pair_to[k]->type_i; base_pair[j].num++; } } index = new int[n_base]; for (i = 0; i < n_base; i++) index[i] = 0; num_pair = (n_base + 1) * 3; best_pairs = new int[num_pair]; num_pair = 0; for (i = 0; i < n_base; i++) { if (index[i]) continue; pair_type = -1; j = get_best_pair(&(base_pair[i]), pair_type, type_i, type_j, cis_or_trans, index); if (j >= 0) { k = get_best_pair(&(base_pair[j]), pair_type, type_i, type_j, cis_or_trans, index); if (k == i) { index[i] = 1; index[j] = 1; best_pairs[6 * num_pair] = i; best_pairs[6*num_pair+1] = j; best_pairs[6*num_pair+2] = pair_type; best_pairs[6*num_pair+3] = type_i; best_pairs[6*num_pair+4] = type_j; best_pairs[6*num_pair+5] = cis_or_trans; /* if (pair_type == 2) best_pairs[3*num_pair+2] = 1; else best_pairs[3*num_pair+2] = 0; */ num_pair++; } } } if (num_pair) re_ordering(num_pair, best_pairs, _base); delete [] best_pairs; delete [] index; for (i = 0; i < n_base; i++) { if (base_pair[i].pair_to) { for (int j = 0; j < base_pair[i].num; j++) if (base_pair[i].pair_to[j]) delete base_pair[i].pair_to[j]; delete [] base_pair[i].pair_to; } } delete [] base_pair; return num_pair; } void _VALIDATION::insert_base_pairs(const int n_contact, CONTACT **contact) { BTreeObj<_BPIndex> _BPIdx; _BPIndex _bpidx; char symmetry[10], bspair_type[35]; int i = 0, ret = 0; for (i = 0; i < n_contact; i++) { _bpidx.set_a_chnid(contact[i]->a_chnid); _bpidx.set_a_resnum(contact[i]->a_resnum); _bpidx.set_a_resnam(contact[i]->a_resnam); _bpidx.set_a_atomtyp(contact[i]->a_atomtyp); _bpidx.set_b_chnid(contact[i]->b_chnid); _bpidx.set_b_resnum(contact[i]->b_resnum); _bpidx.set_b_resnam(contact[i]->b_resnam); _bpidx.set_b_atomtyp(contact[i]->b_atomtyp); _bpidx.set_pos(i); _BPIdx.Add(_bpidx); } _BPIdx.GoFirst(); while (ret >= 0) { _bpidx = _BPIdx.GetCurrData(); i = _bpidx.pos(); symmetry[0] = '\0'; if (contact[i]->sym > 1) { sprintf(symmetry, "%d_%d%d%d", contact[i]->sym, contact[i]->transl[0]+5, contact[i]->transl[1]+5, contact[i]->transl[2]+5); } if (contact[i]->type > 0 && contact[i]->type < 30) strcpy(bspair_type, _base_pair_type[contact[i]->type]); else if (!res_name_cmp_na(contact[i]->a_resnam, "A") && (!res_name_cmp_na(contact[i]->b_resnam, "T") || !res_name_cmp_na(contact[i]->b_resnam, "U")) || (!res_name_cmp_na(contact[i]->a_resnam, "T") || !res_name_cmp_na(contact[i]->a_resnam, "U")) && !res_name_cmp_na(contact[i]->b_resnam, "A") || !res_name_cmp_na(contact[i]->a_resnam, "G") && !res_name_cmp_na(contact[i]->b_resnam, "C") || !res_name_cmp_na(contact[i]->a_resnam, "C") && !res_name_cmp_na(contact[i]->b_resnam, "G")) sprintf(bspair_type, "%s-%s PAIR", contact[i]->a_resnam, contact[i]->b_resnam); else sprintf(bspair_type, "%s-%s MISPAIR", contact[i]->a_resnam, contact[i]->b_resnam); Insert_a_Bspair(contact[i]->a_atomtyp, contact[i]->a_atomnum, contact[i]->a_resnum, contact[i]->b_atomtyp, contact[i]->b_atomnum, contact[i]->b_resnum, symmetry, bspair_type); // if (strstr(bspair_type, "MISPAIR") /* && contact[i]->sym > 1 */) { if (!(!res_name_cmp_na(contact[i]->a_resnam, "A") && (!res_name_cmp_na(contact[i]->b_resnam, "T") || !res_name_cmp_na(contact[i]->b_resnam, "U")) || (!res_name_cmp_na(contact[i]->a_resnam, "T") || !res_name_cmp_na(contact[i]->a_resnam, "U")) && !res_name_cmp_na(contact[i]->b_resnam, "A") || !res_name_cmp_na(contact[i]->a_resnam, "G") && !res_name_cmp_na(contact[i]->b_resnam, "C") || !res_name_cmp_na(contact[i]->a_resnam, "C") && !res_name_cmp_na(contact[i]->b_resnam, "G"))) { Insert_a_Mismat(contact[i]->a_atomnum, contact[i]->a_resnum, contact[i]->b_atomnum, contact[i]->b_resnum); } ret = _BPIdx.GoToNext(); } } void _VALIDATION::re_ordering(const int num_pair, int *best_pairs, _BASE **_base) { int i = 0, i1 = 0, i2 = 0, ip = 0, j = 0, j1 = 0, j2 = 0, k = 0, m = 0, ib = 0; int n = 0, num_ends = 0, direction[6], ddidx[CNUM]; int m1 = 0, k0 = 0, k2 = 0, k3 = 0, wrong_link = 0, anti_parallel, parallel; int **end_list = NULL, **bp_order = NULL, *matched_idx = NULL; int **helix_idx = NULL, *bp_idx = NULL, *index = NULL; double ddmin[CNUM], temp = -9999.9, di1_i2, di1_j2, dj1_i2, dj1_j2, d, d2, d3; COORD zave, txyz, txyz2, *bp_xyz = NULL; end_list = new int*[num_pair]; bp_order = new int*[num_pair]; bp_xyz = new COORD[num_pair]; for (i = 0; i < num_pair; i++) { end_list[i] = new int[3]; bp_order[i] = new int[3]; vector_middle(bp_xyz[i], _base[best_pairs[6 * i]]->org, _base[best_pairs[6*i+1]]->org); for (j = 0; j < 3; j++) { end_list[i][j] = -1; bp_order[i][j] = -1; } } for (i = 0; i < num_pair; i++) { for (j = 0; j < CNUM; j++) { ddmin[j] = XBIG; ddidx[j] = -1; } di1_i2 = vector_dot_product(_base[best_pairs[6 * i]]->orien[2], _base[best_pairs[6*i+1]]->orien[2]); if (di1_i2 <= 0.0) vector_difference(zave, _base[best_pairs[6 * i]]->orien[2], _base[best_pairs[6*i+1]]->orien[2]); else vector_sum(zave, _base[best_pairs[6 * i]]->orien[2], _base[best_pairs[6*i+1]]->orien[2]); for (j = 0; j < num_pair; j++) { if (i != j) { vector_difference(txyz, bp_xyz[i], bp_xyz[j]); temp = vector_length(txyz); for (k = 0; k < CNUM; k++) { if (temp < ddmin[k]) { for (m = CNUM - 1; m > k; m--) { if (ddidx[m-1] >= 0) { ddidx[m] = ddidx[m-1]; ddmin[m] = ddmin[m-1]; } } ddmin[k] = temp; ddidx[k] = j; break; } } } } if (ddidx[0] >= 0 && ddidx[1] >= 0) { if (ddmin[0] > HELIX_CHG) { end_list[num_ends++][0] = i; n = 2; } else { vector_difference(txyz, bp_xyz[i], bp_xyz[ddidx[0]]); d = vector_dot_product(zave, txyz); if (ddidx[2] >= 0 && ddmin[1] <= HELIX_CHG && ddmin[2] <= HELIX_CHG) { vector_difference(txyz, bp_xyz[i], bp_xyz[ddidx[1]]); vector_difference(txyz2, bp_xyz[i], bp_xyz[ddidx[2]]); d2 = vector_dot_product(zave, txyz); d3 = vector_dot_product(zave, txyz2); if (d * d2 < 0.0 && d * d3 < 0.0 && fabs(d2) > fabs(d3)) { k = ddidx[1]; ddidx[1] = ddidx[2]; ddidx[2] = k; temp = ddmin[1]; ddmin[1] = ddmin[2]; ddmin[2] = temp; } } n = -1; for (j = 1; j < CNUM && ddidx[j] >= 0; j++) { if (ddmin[j] > HELIX_CHG) break; vector_difference(txyz2, bp_xyz[i], bp_xyz[ddidx[j]]); d2 = vector_dot_product(zave, txyz2); if (d * d2 < 0.0) { n = j; bp_order[i][0] = -2; bp_order[i][1] = ddidx[0]; bp_order[i][2] = ddidx[j]; break; } } if (n < 0) { n = 2; end_list[num_ends][0] = i; end_list[num_ends][1] = ddidx[0]; bp_order[i][1] = ddidx[0]; vector_difference(txyz2, bp_xyz[ddidx[1]], bp_xyz[ddidx[0]]); d2 = vector_dot_product(zave, txyz2); if (d * d2 < 0.0 && vector_length(txyz2) <= HELIX_CHG) { end_list[num_ends][2] = ddidx[1]; bp_order[i][2] = ddidx[1]; } num_ends++; } } } } if (!num_ends) { end_list[num_ends][0] = 0; if (num_pair == 2) { if (temp <= HELIX_CHG) { end_list[num_ends][1] = 1; num_ends++; end_list[num_ends][0] = 1; end_list[num_ends][1] = 0; num_ends++; } else { num_ends++; end_list[num_ends][0] = 1; num_ends++; } } else num_ends++; } matched_idx = new int[num_pair]; bp_idx = new int[num_pair]; for (i = 0; i < num_pair; i++) { matched_idx[i] = 0; bp_idx[i] = i; } helix_idx = new int*[num_ends]; for (i = 0; i < num_ends; i++) { helix_idx[i] = new int[2]; helix_idx[i][0] = helix_idx[i][1] = 0; } helix_idx[0][0] = 0; n = 0; for (i = 0; i < num_ends && ip < num_pair; i++) { // if (end_list[i][0] == -1) continue; k = 0; k0 = 0; for (j = 0; j < 3; j++) { if (end_list[i][j] >= 0) { k += matched_idx[end_list[i][j]]; k0++; } } if (k == k0) continue; for (j = 0; j < 3 && ip < num_pair; j++) { if (end_list[i][j] >= 0 && !matched_idx[end_list[i][j]]) { bp_idx[ip] = end_list[i][j]; matched_idx[bp_idx[ip]] = 1; ip++; } } for (j = 0; j < num_pair && ip < num_pair; j++) { if (ip < 1) break; k = bp_idx[ip - 1]; k2 = bp_order[k][1]; k3 = bp_order[k][2]; if (bp_order[k][0] == -1) { if (k2 >= 0 && !matched_idx[k2] && k3 == -1) { bp_idx[ip++] = k2; matched_idx[k2] = 1; } break; } m = matched_idx[k2] + matched_idx[k3]; if (m == 2 || m == 0) break; if (k2 == bp_idx[ip - 2]) { bp_idx[ip++] = k3; matched_idx[k3] = 1; } else if (k3 == bp_idx[ip - 2]) { bp_idx[ip++] = k2; matched_idx[k2] = 1; } else break; } helix_idx[n][1] = ip; if (ip && ip < num_pair) { if ((n+1) >= num_ends) { for (k = 0; k < num_pair; k++) { delete [] end_list[k]; delete [] bp_order[k]; } delete [] end_list; delete [] bp_order; delete [] bp_xyz; delete [] matched_idx; delete [] bp_idx; for (k = 0; k < num_ends; k++) delete [] helix_idx[k]; delete [] helix_idx; return; } helix_idx[++n][0] = ip; } } if (helix_idx[n][1] == 0) n--; wrong_link = 0; for (i = 0; i <= n; i++) { for (j = helix_idx[i][0]; j < helix_idx[i][1]; j++) { ib = best_pairs[6 * bp_idx[j]]; if (_base[ib]->exist_P && _base[ib]->exist_O) { di1_i2 = distance_ab(_base[ib], 1, _base[ib], 0); if (di1_i2 <= O3P_UPPER) { wrong_link = 1; break; } } } if (wrong_link) break; } if (!wrong_link) { index = new int[num_pair]; for (i = 0; i < num_pair; i++) index[i] = 0; for (i = 0; i <= n; i++) { if ((helix_idx[i][1] - 1) > helix_idx[i][0]) { j = helix_idx[i][0]; m = bp_idx[j]; m1 = bp_idx[j+1]; get_ij(m, index[m], best_pairs, i1, j1); get_ij(m1, index[m1], best_pairs, i2, j2); k = is_linked(_base[i1], _base[i2]); if (k == -1) index[m] = !index[m]; else if (!k) { reverse(helix_idx[i][0], helix_idx[i][1] - helix_idx[i][0], bp_idx); m = bp_idx[j]; m1 = bp_idx[j+1]; get_ij(m, index[m], best_pairs, i1, j1); get_ij(m1, index[m1], best_pairs, i2, j2); k = is_linked(_base[i1], _base[i2]); if (k == -1) index[m] = !index[m]; else if (!k) reverse(helix_idx[i][0], helix_idx[i][1] - helix_idx[i][0], bp_idx); } } for (j = helix_idx[i][0]; j < helix_idx[i][1] - 1; j++) { m = bp_idx[j]; m1 = bp_idx[j+1]; get_ij(m, index[m], best_pairs, i1, j1); get_ij(m1, index[m1], best_pairs, i2, j2); k = 0; if (is_wc_geometry(_base[i1], _base[j1]) && is_wc_geometry(_base[i2], _base[j2])) { vector_sum(txyz, _base[i1]->orien[0], _base[j1]->orien[0]); vector_sum(txyz2, _base[i2]->orien[0], _base[j2]->orien[0]); if (magang(txyz, txyz2) > 125.0 || is_linked(_base[i1], _base[i2]) || is_linked(_base[j1], _base[j2])) k = 0; else { vector_difference(txyz, _base[i1]->orien[2], _base[j1]->orien[2]); vector_difference(txyz2, _base[i2]->orien[2], _base[j2]->orien[2]); vector_normalize(txyz); vector_normalize(txyz2); di1_i2 = vector_dot_product(txyz, txyz2); vector_difference(txyz2, _base[j2]->orien[2], _base[i2]->orien[2]); vector_normalize(txyz2); di1_j2 = vector_dot_product(txyz, txyz2); if (di1_i2 < 0.0 && di1_j2 > 0.0) k = 1; } } if (k) index[m1] = !index[m1]; else { if (_base[i1]->exist_O && _base[j1]->exist_O && _base[i2]->exist_O && _base[j2]->exist_O) { di1_i2 = distance_ab(_base[i1], 0, _base[i2], 0); di1_j2 = distance_ab(_base[i1], 0, _base[j2], 0); dj1_i2 = distance_ab(_base[j1], 0, _base[i2], 0); dj1_j2 = distance_ab(_base[j1], 0, _base[j2], 0); if (di1_i2 > di1_j2 && dj1_j2 > dj1_i2) k = 1; } if (k) index[m1] = !index[m1]; else { if (!is_linked(_base[i1], _base[i2]) && !is_linked(_base[j1], _base[j2]) && (is_linked(_base[i1], _base[j2]) || is_linked(_base[j1], _base[i2]))) k = 1; if (k) index[m1] = !index[m1]; else { k = check_others(_base, i1, j1, i2, j2); if (k) index[m1] = !index[m1]; } } } get_ij(m, index[m], best_pairs, i1, j1); get_ij(m1, index[m1], best_pairs, i2, j2); if (is_linked(_base[i1], _base[i2]) == -1) index[m1] = !index[m1]; } if (check_direction(helix_idx[i], bp_idx, index, best_pairs, _base, direction)) { if (direction[0] + direction[1] + direction[3] + direction[4]) { for (j = helix_idx[i][0]; j < helix_idx[i][1] - 1; j++) { m = bp_idx[j]; get_ij(m, index[m], best_pairs, i1, j1); m1 = bp_idx[j+1]; get_ij(m1, index[m1], best_pairs, i2, j2); if (!is_linked(_base[i1], _base[i2]) && !is_linked(_base[j1], _base[j2]) && ((is_linked(_base[i1], _base[j2]) == 1) || (is_linked(_base[i1], _base[j2]) && is_linked(_base[j1], _base[i2])))) index[m1] = !index[m1]; } } } else { anti_parallel = direction[0] > direction[1] && direction[3] < direction[4]; parallel = direction[0] > direction[1] && direction[3] > direction[4]; for (j = helix_idx[i][0]; j < helix_idx[i][1] - 1; j++) { m = bp_idx[j]; get_ij(m, index[m], best_pairs, i1, j1); m1 = bp_idx[j+1]; get_ij(m1, index[m1], best_pairs, i2, j2); k = is_linked(_base[j1], _base[j2]); if (!is_linked(_base[i1], _base[i2]) && ((anti_parallel && k == 1) || (parallel && (k == -1)))) index[m1] = !index[m1]; get_ij(m1, index[m1], best_pairs, i2, j2); if (!is_linked(_base[i1], _base[i2]) && !is_linked(_base[j1], _base[j2])) { if ((anti_parallel && is_linked(_base[j1], _base[i2]) == 1) || (parallel && is_linked(_base[i1], _base[j2]) == -1)) index[m] = !index[m]; else if ((parallel && is_linked(_base[j1], _base[i2]) == -1) || (anti_parallel && is_linked(_base[i1], _base[j2]) == 1)) index[m1] = !index[m1]; } } } check_direction(helix_idx[i], bp_idx, index, best_pairs, _base, direction); for (j = helix_idx[i][0]; j < helix_idx[i][1]; j++) { m = bp_idx[j]; if (index[m]) { i1 = best_pairs[6 * m]; best_pairs[6 * m] = best_pairs[6*m+1]; best_pairs[6*m+1] = i1; i1 = best_pairs[6*m+3]; best_pairs[6*m+3] = best_pairs[6*m+4]; best_pairs[6*m+4] = i1; } } } delete [] index; } for (i = 0; i <= n; i++) { i1 = helix_idx[i][1] - helix_idx[i][0]; index = new int[i1 * 6]; for (i2 = 0; i2 < i1; i2++) { j = i2 + helix_idx[i][0]; j1 = bp_idx[j]; index[6 * i2] = best_pairs[6 * j1]; index[6*i2+1] = best_pairs[6*j1+1]; index[6*i2+2] = best_pairs[6*j1+2]; index[6*i2+3] = best_pairs[6*j1+3]; index[6*i2+4] = best_pairs[6*j1+4]; index[6*i2+5] = best_pairs[6*j1+5]; } get_parameters(i1, index, _base); delete [] index; } for (i = 0; i < num_pair; i++) { delete [] end_list[i]; delete [] bp_order[i]; } delete [] end_list; delete [] bp_order; delete [] bp_xyz; delete [] matched_idx; delete [] bp_idx; for (i = 0; i < num_ends; i++) delete [] helix_idx[i]; delete [] helix_idx; } void _VALIDATION::get_parameters(const int num_bp, int *pair_num, _BASE **_base) { int ds = 2, i, j, xdir, helix_value; double dsum, spi[6], hpi[6], *twist = NULL; COORD mfi[3], mfoi, mfi_prev[3], mfoi_prev, *aveS = NULL, *aveH = NULL; COORD mst_orien[3], mst_org, mst_orienH[3], mst_orgH; char sym[20]; _RESIDUE *res1 = NULL, *res2 = NULL, *res3 = NULL, *res4 = NULL; for (j = 0; j < num_bp; j++) { /* new change */ if (vector_dot_product(_base[pair_num[6 * j]]->orien[2], _base[pair_num[6 * j + 1]]->orien[2]) < 0.0) { vector_reverse(_base[pair_num[6 * j + 1]]->orien[1]); vector_reverse(_base[pair_num[6 * j + 1]]->orien[2]); } } xdir = 0; for (i = 0; i < ds; i++) { for (j = 0; j < num_bp; j++) { if (j < (num_bp - 1)) vector_difference(mst_org, _base[pair_num[6 * (j+1) + i]]->org, _base[pair_num[6 * j + i]]->org); if (vector_dot_product(mst_org, _base[pair_num[6 * j + i]]->orien[2]) < 0.0) xdir++; } } if (xdir == ds * num_bp) { for (i = 0; i < ds; i++) { for (j = 0; j < num_bp; j++) { vector_reverse(_base[pair_num[6 * j + i]]->orien[0]); vector_reverse(_base[pair_num[6 * j + i]]->orien[2]); } } } dsum = 0; twist = new double[num_bp - 1]; aveS = new COORD[num_bp - 1]; aveH = new COORD[num_bp - 1]; for (i = 0; i < num_bp - 1; i++) { vector_init(aveS[i], XBIG, XBIG, XBIG); vector_init(aveH[i], XBIG, XBIG, XBIG); } for (i = 0; i < num_bp; i++) { bpstep_par(_base[pair_num[6*i+1]]->orien, _base[pair_num[6*i+1]]->org, _base[pair_num[6 * i]]->orien, _base[pair_num[6 * i]]->org, spi, mfi, mfoi); res3 = _base[pair_num[6 * i]]->res; res4 = _base[pair_num[6*i+1]]->res; if (_n_bb_params >= _alloc_bb_params) { _alloc_bb_params += 20; BASE_BASE_PARAMS **tmp = _bb_params; _bb_params = new BASE_BASE_PARAMS*[_alloc_bb_params]; for (j = 0; j < _n_bb_params; j++) _bb_params[j] = tmp[j]; for (j = _n_bb_params; j < _alloc_bb_params; j++) _bb_params[j] = NULL; delete [] tmp; } _bb_params[_n_bb_params] = new BASE_BASE_PARAMS; _bb_params[_n_bb_params]->ChnID_I.Copy(res3->GetFirstAtom()->chnid()); strcpy(_bb_params[_n_bb_params]->ChnID_I_PDB, res3->GetFirstAtom()->pdb_chnid()); strcpy(_bb_params[_n_bb_params]->Base_I, res3->ResName()); strcpy(_bb_params[_n_bb_params]->Sym_I, "1_555"); _bb_params[_n_bb_params]->ResNum_I = res3->res_no(); _bb_params[_n_bb_params]->ResNum_I_PDB = res3->pdb_res_no(); strcpy(_bb_params[_n_bb_params]->InsCode_I, res3->GetFirstAtom()->ins_code()); _bb_params[_n_bb_params]->ChnID_J.Copy(res4->GetFirstAtom()->chnid()); strcpy(_bb_params[_n_bb_params]->ChnID_J_PDB, res4->GetFirstAtom()->pdb_chnid()); strcpy(_bb_params[_n_bb_params]->Base_J, res4->ResName()); if (_base[pair_num[6 * i + 1]]->op || _base[pair_num[6 * i + 1]]->lx || _base[pair_num[6 * i + 1]]->ly || _base[pair_num[6 * i + 1]]->lz) { sprintf(sym, "%d_%d%d%d", _base[pair_num[6 * i + 1]]->op + 1, _base[pair_num[6 * i + 1]]->lx + 5, _base[pair_num[6 * i + 1]]->ly + 5, _base[pair_num[6 * i + 1]]->lz + 5); strcpy(_bb_params[_n_bb_params]->Sym_J, sym); } else strcpy(_bb_params[_n_bb_params]->Sym_J, "1_555"); _bb_params[_n_bb_params]->ResNum_J = res4->res_no(); _bb_params[_n_bb_params]->ResNum_J_PDB = res4->pdb_res_no(); strcpy(_bb_params[_n_bb_params]->InsCode_J, res4->GetFirstAtom()->ins_code()); _bb_params[_n_bb_params]->Shear = spi[0]; _bb_params[_n_bb_params]->Stretch = spi[1]; _bb_params[_n_bb_params]->Stagger = spi[2]; _bb_params[_n_bb_params]->Buckle = spi[3]; _bb_params[_n_bb_params]->Propel = spi[4]; _bb_params[_n_bb_params]->Opening = spi[5]; if (pair_num[6*i+2] > 0 && pair_num[6*i+2] < 30) /* strcpy(_bb_params[_n_bb_params]->_saenger_classification, _saenger_classification[pair_num[6 * i + 2]]); */ sprintf(_bb_params[_n_bb_params]->_saenger_classification, "%d", pair_num[6*i+2]); else strcpy(_bb_params[_n_bb_params]->_saenger_classification, ""); if (pair_num[6*i+3] && pair_num[6*i+4] && pair_num[6*i+5]) { /* sprintf(_bb_params[_n_bb_params]->_leontis_westhof_classification, "%s %s/%s", _gb_dir[pair_num[6*i+5]], _edges[pair_num[6*i+3]], _edges[pair_num[6*i+4]]); */ sprintf(_bb_params[_n_bb_params]->_leontis_westhof_classification, "%d", _leontis_westhof_index[pair_num[6*i+5]-1][pair_num[6*i+3]-1][pair_num[6*i+4]-1]); } else if (pair_num[6*i+2] == HBOND_TYPE_19 || pair_num[6*i+2] == HBOND_TYPE_20) /* strcpy(_bb_params[_n_bb_params]->_leontis_westhof_classification, "Cis Watson-Crick/Watson-Crick"); */ strcpy(_bb_params[_n_bb_params]->_leontis_westhof_classification, "1"); else strcpy(_bb_params[_n_bb_params]->_leontis_westhof_classification, ""); _n_bb_params++; if (!i) { mfoi_prev = mfoi; for (j = 0; j < 3; j++) mfi_prev[j] = mfi[j]; continue; } bpstep_par(mfi_prev, mfoi_prev, mfi, mfoi, spi, mst_orien, mst_org); helical_par(mfi_prev, mfoi_prev, mfi, mfoi, hpi, mst_orienH, mst_orgH); mfoi_prev = mfoi; for (j = 0; j < 3; j++) mfi_prev[j] = mfi[j]; if (_base[pair_num[6 * i]]->exist_P && _base[pair_num[6*(i-1)+1]]->exist_P) { project_xyzP(_base[pair_num[6 * i]]->P_xyz, _base[pair_num[6*(i-1)+1]]->P_xyz, mst_orien, mst_org, aveS[i-1]); project_xyzP(_base[pair_num[6 * i]]->P_xyz, _base[pair_num[6*(i-1)+1]]->P_xyz, mst_orienH, mst_orgH, aveH[i-1]); } dsum += spi[5]; twist[i - 1] = spi[5]; res1 = _base[pair_num[6 * (i-1)]]->res; res2 = _base[pair_num[6*(i-1)+1]]->res; if (_n_ib_params >= _alloc_ib_params) { _alloc_ib_params += 20; INTER_BASE_PARAMS **tmp = _ib_params; _ib_params = new INTER_BASE_PARAMS*[_alloc_ib_params]; for (j = 0; j < _n_ib_params; j++) _ib_params[j] = tmp[j]; for (j = _n_ib_params; j < _alloc_ib_params; j++) _ib_params[j] = NULL; delete [] tmp; } _ib_params[_n_ib_params] = new INTER_BASE_PARAMS; _ib_params[_n_ib_params]->ChnID_I1.Copy(res1->GetFirstAtom()->chnid()); strcpy(_ib_params[_n_ib_params]->ChnID_I1_PDB, res1->GetFirstAtom()->pdb_chnid()); strcpy(_ib_params[_n_ib_params]->Base_I1, res1->ResName()); strcpy(_ib_params[_n_ib_params]->Sym_I1, "1_555"); _ib_params[_n_ib_params]->ResNum_I1 = res1->res_no(); _ib_params[_n_ib_params]->ResNum_I1_PDB = res1->pdb_res_no(); strcpy(_ib_params[_n_ib_params]->InsCode_I1, res1->GetFirstAtom()->ins_code()); _ib_params[_n_ib_params]->ChnID_J1.Copy(res2->GetFirstAtom()->chnid()); strcpy(_ib_params[_n_ib_params]->ChnID_J1_PDB, res2->GetFirstAtom()->pdb_chnid()); strcpy(_ib_params[_n_ib_params]->Base_J1, res2->ResName()); if (_base[pair_num[6*(i-1)+1]]->op || _base[pair_num[6*(i-1)+1]]->lx || _base[pair_num[6*(i-1)+1]]->ly || _base[pair_num[6*(i-1)+1]]->lz) { sprintf(sym, "%d_%d%d%d", _base[pair_num[6*(i-1)+1]]->op + 1, _base[pair_num[6*(i-1)+1]]->lx + 5, _base[pair_num[6*(i-1)+1]]->ly + 5, _base[pair_num[6*(i-1)+1]]->lz + 5); strcpy(_ib_params[_n_ib_params]->Sym_J1, sym); } else strcpy(_ib_params[_n_ib_params]->Sym_J1, "1_555"); _ib_params[_n_ib_params]->ResNum_J1 = res2->res_no(); _ib_params[_n_ib_params]->ResNum_J1_PDB = res2->pdb_res_no(); strcpy(_ib_params[_n_ib_params]->InsCode_J1, res2->GetFirstAtom()->ins_code()); _ib_params[_n_ib_params]->ChnID_I2.Copy(res3->GetFirstAtom()->chnid()); strcpy(_ib_params[_n_ib_params]->ChnID_I2_PDB, res3->GetFirstAtom()->pdb_chnid()); strcpy(_ib_params[_n_ib_params]->Base_I2, res3->ResName()); strcpy(_ib_params[_n_ib_params]->Sym_I2, "1_555"); _ib_params[_n_ib_params]->ResNum_I2 = res3->res_no(); _ib_params[_n_ib_params]->ResNum_I2_PDB = res3->pdb_res_no(); strcpy(_ib_params[_n_ib_params]->InsCode_I2, res3->GetFirstAtom()->ins_code()); _ib_params[_n_ib_params]->ChnID_J2.Copy(res4->GetFirstAtom()->chnid()); strcpy(_ib_params[_n_ib_params]->ChnID_J2_PDB, res4->GetFirstAtom()->pdb_chnid()); strcpy(_ib_params[_n_ib_params]->Base_J2, res4->ResName()); if (_base[pair_num[6*i+1]]->op || _base[pair_num[6*i+1]]->lx || _base[pair_num[6*i+1]]->ly || _base[pair_num[6*i+1]]->lz) { sprintf(sym, "%d_%d%d%d", _base[pair_num[6*i+1]]->op + 1, _base[pair_num[6*i+1]]->lx + 5, _base[pair_num[6*i+1]]->ly + 5, _base[pair_num[6*i+1]]->lz + 5); strcpy(_ib_params[_n_ib_params]->Sym_J2, sym); } else strcpy(_ib_params[_n_ib_params]->Sym_J2, "1_555"); _ib_params[_n_ib_params]->ResNum_J2 = res4->res_no(); _ib_params[_n_ib_params]->ResNum_J2_PDB = res4->pdb_res_no(); strcpy(_ib_params[_n_ib_params]->InsCode_J2, res4->GetFirstAtom()->ins_code()); _ib_params[_n_ib_params]->Shift = spi[0]; _ib_params[_n_ib_params]->Slide = spi[1]; _ib_params[_n_ib_params]->Rise = spi[2]; _ib_params[_n_ib_params]->Tilt = spi[3]; _ib_params[_n_ib_params]->Roll = spi[4]; _ib_params[_n_ib_params]->Twist = spi[5]; _ib_params[_n_ib_params]->X_disp = hpi[0]; _ib_params[_n_ib_params]->Y_disp = hpi[1]; _ib_params[_n_ib_params]->H_rise = hpi[2]; _ib_params[_n_ib_params]->Incl = hpi[3]; _ib_params[_n_ib_params]->Tip = hpi[4]; _ib_params[_n_ib_params]->H_twist = hpi[5]; _n_ib_params++; } dsum /= (double) (num_bp - 1); if (dsum > 10.0 && dsum < 60.0) { helix_value = DOUBLE_HELIX; for (i = 0; i < num_bp - 1; i++) { if (twist[i] > 0.0 && aveS[i].x > -5.0 && aveS[i].x < -0.5 && aveS[i].y > 7.5 && aveS[i].y < 10.0 && aveS[i].z > -2.0 && aveS[i].z < 3.5 && aveH[i].x > -11.5 && aveH[i].x < 2.5 && aveH[i].y > 1.5 && aveH[i].y < 10.0 && aveH[i].z > -3.0 && aveH[i].z < 9.0) { if (aveS[i].z >= 1.5) { if (helix_value == A_DOUBLE_HELIX) _classification |= A_DOUBLE_HELIX; helix_value = A_DOUBLE_HELIX; } else if (aveS[i].z <= 0.5) { if (helix_value == B_DOUBLE_HELIX) _classification |= B_DOUBLE_HELIX; helix_value = B_DOUBLE_HELIX; } else { _classification |= DOUBLE_HELIX; helix_value = DOUBLE_HELIX; } } else { _classification |= DOUBLE_HELIX; helix_value = DOUBLE_HELIX; } } } else if (dsum < 0) { if (xdir == (ds * num_bp)) _classification |= Z_DOUBLE_HELIX; else _classification |= DOUBLE_HELIX; } else _classification |= DOUBLE_HELIX; delete [] twist; delete [] aveS; delete [] aveH; } void _VALIDATION::create_blocks(_CHAIN *chain, const int color_type) { int Y_or_R, nmatch, i, j, k, l; double eRing_xyz[27], sRing_xyz[27], rt[12], ave_exyz[3], ave_sxyz[3]; double *block = NULL; _BLOCK **tmp = NULL; _RESIDUE *res = chain->GetFirstResidue(); while (res) { Y_or_R = 0; nmatch = find_match(res, eRing_xyz, sRing_xyz, &block, Y_or_R); if (nmatch > 2) { superposition_matrix(nmatch, eRing_xyz, sRing_xyz, rt, ave_exyz, ave_sxyz); for (k = 0; k < 3; k++) { for (l = 0; l < 3; l++) ave_exyz[k] -= rt[k * 3 + l] * ave_sxyz[l]; } if (_n_blocks >= _alloc_blocks) { _alloc_blocks += chain->TotalResidueNumbers(); tmp = _blocks; _blocks = new _BLOCK*[_alloc_blocks]; for (i = 0; i < _n_blocks; i++) _blocks[i] = tmp[i]; for (i = _n_blocks; i < _alloc_blocks; i++) _blocks[i] = NULL; delete [] tmp; } _blocks[_n_blocks] = new _BLOCK; for (j = 0; j < 8; j++) { for (k = 0; k < 3; k++) { ave_sxyz[k] = ave_exyz[k]; for (int l = 0; l < 3; l++) ave_sxyz[k] += rt[k * 3 + l] * block[j * 3 + l]; } _blocks[_n_blocks]->coord[j].x = ave_sxyz[0]; _blocks[_n_blocks]->coord[j].y = ave_sxyz[1]; _blocks[_n_blocks]->coord[j].z = ave_sxyz[2]; } _blocks[_n_blocks]->index = ndb_get_color_index(color_type, res->ResName(), chain->chain_type(), chain->icolor()); _n_blocks++; } res = chain->GetNextResidue(); } } static void realloc_base(_BASE*** _base, const int n_base, int &alloc_base) { int i = 0; _BASE **tmp = NULL; if (n_base >= alloc_base) { alloc_base += 20; tmp = *_base; *_base = new _BASE*[alloc_base]; for (i = 0; i < n_base; i++) (*_base)[i] = tmp[i]; for (i = n_base; i < alloc_base; i++) (*_base)[i] = NULL; delete [] tmp; } } static _BASE *base_frame(_RESIDUE *res, const int op, const int lx, const int ly, const int lz, const NDBCELLCONST *cell) { int k = 0, l = 0, Y_or_R = 0; double *block = NULL, eRing_xyz[27], sRing_xyz[27], rt[12]; double ave_exyz[4], ave_sxyz[4]; int nmatch = find_match(res, eRing_xyz, sRing_xyz, &block, Y_or_R); if (nmatch < 3) return NULL; _BASE *base = new _BASE; base->res = res; _ATOM *atom = res->find_atom("P"); if (atom) { base->exist_P = 1; base->P_xyz = atom->orig(); } else base->exist_P = 0; atom = res->find_atom("O3*"); if (atom) { base->exist_O = 1; base->O_xyz = atom->orig(); } else base->exist_O = 0; if (Y_or_R) atom = res->find_atom("N9"); else atom = res->find_atom("N1"); if (atom) { base->exist_N = 1; base->N_xyz = atom->orig(); } else base->exist_N = 0; atom = res->find_atom("C1*"); if (atom) { base->exist_C = 1; base->C_xyz = atom->orig(); } else base->exist_C = 0; base->op = op; base->lx = lx; base->ly = ly; base->lz = lz; if (op || lx || ly || lz) { for (k = 0; k < nmatch; k++) { symmetry_operation(&eRing_xyz[3*k], &eRing_xyz[3*k+1], &eRing_xyz[3*k+2], cell, op, lx, ly, lz); } if (base->exist_P) symmetry_operation(base->P_xyz, cell, op, lx, ly, lz); if (base->exist_O) symmetry_operation(base->O_xyz, cell, op, lx, ly, lz); if (base->exist_N) symmetry_operation(base->N_xyz, cell, op, lx, ly, lz); if (base->exist_C) symmetry_operation(base->C_xyz, cell, op, lx, ly, lz); } superposition_matrix(nmatch, eRing_xyz, sRing_xyz, rt, ave_exyz, ave_sxyz); for (k = 0; k < 3; k++) { for (l = 0; l < 3; l++) ave_exyz[k] -= rt[k * 3 + l] * ave_sxyz[l]; } base->org.x = ave_exyz[0]; base->org.y = ave_exyz[1]; base->org.z = ave_exyz[2]; for (k = 0; k < 3; k++) { base->orien[k].x = rt[k]; base->orien[k].y = rt[3 + k]; base->orien[k].z = rt[6 + k]; } return base; } static int find_match(_RESIDUE *residue, double *eRing_xyz, double *sRing_xyz, double **block, int& Y_or_R) { REF_NA *ref = NULL; const char **RingAtom = NULL; int RingAtom_num = 0, k, l, nmatch, index; char resname[5]; _ATOM *atom = NULL; strcpy(resname, residue->ResName()); get_parent_name(resname); Y_or_R = 0; if (!strcmp(resname, "A")) { ref = &STAND_A; RingAtom = R_RingAtom; RingAtom_num = R_RINGATOM; *block = &R_BLOCK[0][0]; Y_or_R = 1; } else if (!strcmp(resname, "C")) { ref = &STAND_C; RingAtom = Y_RingAtom; RingAtom_num = Y_RINGATOM; *block = &Y_BLOCK[0][0]; } else if (!strcmp(resname, "G")) { ref = &STAND_G; RingAtom = R_RingAtom; RingAtom_num = R_RINGATOM; *block = &R_BLOCK[0][0]; Y_or_R = 1; } else if (!strcmp(resname, "P")) { ref = &STAND_P; RingAtom = Y_RingAtom; RingAtom_num = Y_RINGATOM; *block = &Y_BLOCK[0][0]; } else if (!strcmp(resname, "T")) { ref = &STAND_T; RingAtom = Y_RingAtom; RingAtom_num = Y_RINGATOM; *block = &Y_BLOCK[0][0]; } else if (!strcmp(resname, "U")) { ref = &STAND_U; RingAtom = Y_RingAtom; RingAtom_num = Y_RINGATOM; *block = &Y_BLOCK[0][0]; } if (!ref || !RingAtom) return 0; nmatch = 0; for (k = 0; k < RingAtom_num; k++) { index = find_atom_index(ref->num_atom, ref->AtomName, RingAtom[k]); atom = residue->find_atom(RingAtom[k]); if (index >= 0 && atom) { eRing_xyz[nmatch * 3] = atom->orig().x; eRing_xyz[nmatch*3+1] = atom->orig().y; eRing_xyz[nmatch*3+2] = atom->orig().z; for (l = 0; l < 3; l++) sRing_xyz[nmatch * 3 + l] = ref->ref[index][l]; nmatch++; } } return nmatch; } static _BASE_PAIR *find_all_pairs(const int n_base, _BASE **base) { int i = 0, j = 0, k = 0, i1 = 0, j1 = 0; double d, dd, dn, dv, dp; COORD coord, coord1; _PAIR_TO **tmp = NULL; _BASE_PAIR *base_pair = new _BASE_PAIR[n_base]; for (i = 0; i < n_base; i++) { base_pair[i].num = 0; base_pair[i].alloc_num = 0; base_pair[i].pair_to = NULL; } for (i = 0; i < n_base - 1; i++) { for (j = i + 1; j < n_base; j++) { dd = vector_dot_product(base[i]->orien[2], base[j]->orien[2]); dv = 90.0 - fabs(dot2ang(dd) - 90.0); dn = MIN_N91 + 1.0; if (base[i]->exist_N && base[j]->exist_N) { vector_difference(coord, base[j]->N_xyz, base[i]->N_xyz); dn = vector_length(coord); } vector_difference(coord, base[j]->org, base[i]->org); d = vector_length(coord); if (dd <= 0.0) vector_difference(coord1, base[i]->orien[2], base[j]->orien[2]); else vector_sum(coord1, base[i]->orien[2], base[j]->orien[2]); vector_normalize(coord1); dp = fabs(vector_dot_product(coord, coord1)); if (d <= MAX_DORG && dp <= MAX_DPRJ && dv <= MAX_ANGLE && dn >= MIN_N91) { if (!base[i]->op && !base[i]->lx && !base[i]->ly && !base[i]->lz) { i1 = i; j1 = j; } else { i1 = j; j1 = i; } if (base_pair[i1].num >= base_pair[i1].alloc_num) { base_pair[i1].alloc_num += 2; tmp = base_pair[i1].pair_to; base_pair[i1].pair_to = new _PAIR_TO*[base_pair[i1].alloc_num]; for (k = 0; k < base_pair[i1].num; k++) base_pair[i1].pair_to[k] = tmp[k]; for (k = base_pair[i1].num; k < base_pair[i1].alloc_num; k++) base_pair[i1].pair_to[k] = NULL; delete [] tmp; } base_pair[i1].pair_to[base_pair[i1].num] = new _PAIR_TO; base_pair[i1].pair_to[base_pair[i1].num]->index = j1; base_pair[i1].pair_to[base_pair[i1].num]->d2dp = d + 2.0 * dp; base_pair[i1].pair_to[base_pair[i1].num]->type = -1; base_pair[i1].num++; } } } return base_pair; } static CONTACT **find_real_base_pairs(const int n_base, _BASE **base, _BASE_PAIR *base_pair, int &n_contact, const NDBCELLCONST *cell, _MOLECULE *mol) { int count = 0, i = 0, j = 0, k = 0, l = 0, m = 0, num = 0, reverse = 0, n_alloc = 0; int FstIndex = -1, FstRes = -1, SndIndex = -1, SndRes = -1; int atom_pair_type[100], pos[3], watson_i, hoogsteen_i, sugar_i, watson_j; int hoogsteen_j, sugar_j, max1, max2, max; double dist = 0; char a_resnam[5], b_resnam[5], buffer[15]; _ATOM *a_atom = NULL, *b_atom = NULL; _Atom_Edge ae_tmp, *ext_a = NULL, *ext_b = NULL; ATOM_PAIR atom_pair[100]; BASE_HBOND_TYPE base_tmp, *base_exist = NULL; BTreeObj<_SIndex> _AtomPairIndex; CONTACT **contact = NULL, **c_tmp = NULL; COORD coord, coord1, coord2, cross_coord1, cross_coord2; _DONOR_ACCEPTOR da_tmp, *exist_a = NULL, *exist_b = NULL; _SIndex index; n_contact = 0; for (i = 0; i < n_base - 1; i++) { FstIndex = -1; FstRes = -1; if (!mol->find_ndb_residue(base[i]->res->GetFirstAtom()->chnid(), base[i]->res->ResName(), base[i]->res->GetFirstAtom()->resnum(), "", FstIndex, FstRes)) continue; if (base[i]->res->ResName()[0] == '+') strcpy(a_resnam, &(base[i]->res->ResName()[1])); else strcpy(a_resnam, base[i]->res->ResName()); get_parent_name(a_resnam); if (!strcmp(a_resnam, "P")) strcpy(a_resnam, "U"); strcpy(da_tmp.type, a_resnam); exist_a = (_DONOR_ACCEPTOR *) bsearch(&da_tmp, _base_type, NUM_BASE_TYPE, sizeof(_DONOR_ACCEPTOR), base_type_compare); if (!exist_a) continue; for (k = 0; k < base_pair[i].num; k++) { j = base_pair[i].pair_to[k]->index; SndIndex = -1; SndRes = -1; if (!mol->find_ndb_residue(base[j]->res->GetFirstAtom()->chnid(), base[j]->res->ResName(), base[j]->res->GetFirstAtom()->resnum(), "", SndIndex, SndRes)) continue; if (base[j]->res->ResName()[0] == '+') strcpy(b_resnam, &(base[j]->res->ResName()[1])); else strcpy(b_resnam, base[j]->res->ResName()); get_parent_name(b_resnam); if (!strcmp(b_resnam, "P")) strcpy(b_resnam, "U"); strcpy(da_tmp.type, b_resnam); exist_b = (_DONOR_ACCEPTOR *) bsearch(&da_tmp, _base_type, NUM_BASE_TYPE, sizeof(_DONOR_ACCEPTOR), base_type_compare); if (!exist_b) continue; count = 0; for (l = 0; l < exist_a->da->n_donor; l++) { a_atom = base[i]->res->find_atom(exist_a->da->donor[l]); if (!a_atom) continue; for (m = 0; m < exist_b->da->n_acceptor; m++) { b_atom = base[j]->res->find_atom(exist_b->da->acceptor[m]); if (!b_atom) continue; dist = cal_dist(a_atom, b_atom, cell, base[j]->op, base[j]->lx, base[j]->ly, base[j]->lz); if (dist >= MIN_CUTOFF && dist <= MAX_CUTOFF) { strcpy(atom_pair[count].a_atomtyp, exist_a->da->donor[l]); strcpy(atom_pair[count].b_atomtyp, exist_b->da->acceptor[m]); if ((exist_a->da->donor_type[l] & 2) & (exist_b->da->acceptor_type[m] & 2)) atom_pair_type[count] = 2; else if ((exist_a->da->donor_type[l] & 1) & (exist_b->da->acceptor_type[m] & 1)) atom_pair_type[count] = 1; else atom_pair_type[count] = 0; count++; } } } for (l = 0; l < exist_a->da->n_acceptor; l++) { a_atom = base[i]->res->find_atom(exist_a->da->acceptor[l]); if (!a_atom) continue; for (m = 0; m < exist_b->da->n_donor; m++) { if (!strcmp(exist_a->da->acceptor[l], "O2*") && !strcmp(exist_b->da->donor[m], "O2*")) continue; b_atom = base[j]->res->find_atom(exist_b->da->donor[m]); if (!b_atom) continue; dist = cal_dist(a_atom, b_atom, cell, base[j]->op, base[j]->lx, base[j]->ly, base[j]->lz); if (dist >= MIN_CUTOFF && dist <= MAX_CUTOFF) { strcpy(atom_pair[count].a_atomtyp, exist_a->da->acceptor[l]); strcpy(atom_pair[count].b_atomtyp, exist_b->da->donor[m]); if ((exist_a->da->acceptor_type[l] & 2) & (exist_b->da->donor_type[m] & 2)) atom_pair_type[count] = 2; else if ((exist_a->da->acceptor_type[l] & 1) & (exist_b->da->donor_type[m] & 1)) atom_pair_type[count] = 1; else atom_pair_type[count] = 0; count++; } } } if (!count) continue; for (l = 0; l < count; l++) if (atom_pair_type[l] == 2) break; if (l == count) continue; if (strcmp(a_resnam, b_resnam) > 0) { sprintf(base_tmp.type, "%s%s", b_resnam, a_resnam); reverse = 1; } else { sprintf(base_tmp.type, "%s%s", a_resnam, b_resnam); reverse = 0; } base_exist = (BASE_HBOND_TYPE *) bsearch(&base_tmp, base_hbond_type, NUM_BASE_HBOND_TYPE, sizeof(BASE_HBOND_TYPE), base_hbond_type_compare); if (!base_exist) continue; _AtomPairIndex.Clear(); for (l = 0; l < count; l++) { if (reverse) sprintf(buffer, "%s_%s", atom_pair[l].b_atomtyp, atom_pair[l].a_atomtyp); else sprintf(buffer, "%s_%s", atom_pair[l].a_atomtyp, atom_pair[l].b_atomtyp); index.set_value(buffer); index.set_pos(l); _AtomPairIndex.Add(index); } base_pair[i].pair_to[k]->type = 0; for (l = 0; l < base_exist->btype->num_base_type; l++) { for (m = 0; m < base_exist->btype->base_type[l].num; m++) { sprintf(buffer, "%s_%s", base_exist->btype->base_type[l].atom_pair[m].a_atomtyp, base_exist->btype->base_type[l].atom_pair[m].b_atomtyp); index.set_value(buffer); if (_AtomPairIndex.Seek(index) >= 0) { index = _AtomPairIndex.GetCurrData(); pos[m] = index.pos(); } else break; } if (m == base_exist->btype->base_type[l].num) { base_pair[i].pair_to[k]->type = base_exist->btype->base_type[l].type; num = m; break; } } if (!base_pair[i].pair_to[k]->type) { num = 0; for (l = 0; l < count; l++) { if (atom_pair_type[l] == 2) { pos[num] = l; num++; break; } } } base_pair[i].pair_to[k]->type_i = 0; base_pair[i].pair_to[k]->type_j = 0; base_pair[i].pair_to[k]->cis_or_trans = 0; if (base_pair[i].pair_to[k]->type != HBOND_TYPE_19 && base_pair[i].pair_to[k]->type != HBOND_TYPE_20) { watson_i = 0; hoogsteen_i = 0; sugar_i = 0; watson_j = 0; hoogsteen_j = 0; sugar_j = 0; for (l = 0; l < count; l++) { strcpy(ae_tmp.atom, atom_pair[l].a_atomtyp); ext_a = (_Atom_Edge *) bsearch(&ae_tmp, exist_a->be->_atoms, exist_a->be->num, sizeof(_Atom_Edge), atom_edge_compare); strcpy(ae_tmp.atom, atom_pair[l].b_atomtyp); ext_b = (_Atom_Edge *) bsearch(&ae_tmp, exist_b->be->_atoms, exist_b->be->num, sizeof(_Atom_Edge), atom_edge_compare); if (!ext_a || !ext_b) continue; if (ext_a->edge & 1) watson_i++; if (ext_a->edge & 2) hoogsteen_i++; if (ext_a->edge & 4) sugar_i++; if (ext_b->edge & 1) watson_j++; if (ext_b->edge & 2) hoogsteen_j++; if (ext_b->edge & 4) sugar_j++; } if ((watson_i > 1 || hoogsteen_i > 1 || sugar_i > 1) && (watson_j > 1 || hoogsteen_j > 1 || sugar_j > 1)) { max1 = (watson_i >= hoogsteen_i) ? watson_i : hoogsteen_i; max2 = (hoogsteen_i >= sugar_i) ? hoogsteen_i : sugar_i; max = (max1 >= max2) ? max1 : max2; if (max == watson_i) base_pair[i].pair_to[k]->type_i = WATSON; else if (max == hoogsteen_i) base_pair[i].pair_to[k]->type_i = HOOGSTEEN; else if (max == sugar_i) base_pair[i].pair_to[k]->type_i = SUGAR; max1 = (watson_j >= hoogsteen_j) ? watson_j : hoogsteen_j; max2 = (hoogsteen_j >= sugar_j) ? hoogsteen_j : sugar_j; max = (max1 >= max2) ? max1 : max2; if (max == watson_j) base_pair[i].pair_to[k]->type_j = WATSON; else if (max == hoogsteen_j) base_pair[i].pair_to[k]->type_j = HOOGSTEEN; else if (max == sugar_j) base_pair[i].pair_to[k]->type_j = SUGAR; } if (base[i]->exist_N && base[i]->exist_C && base[j]->exist_N && base[j]->exist_C) { vector_difference(coord1, base[i]->C_xyz, base[i]->N_xyz); vector_difference(coord2, base[j]->C_xyz, base[j]->N_xyz); vector_difference(coord, base[j]->N_xyz, base[i]->N_xyz); vector_cross_product(cross_coord1, coord1, coord); vector_cross_product(cross_coord2, coord2, coord); if (vector_dot_product(cross_coord1, cross_coord2) > 0) base_pair[i].pair_to[k]->cis_or_trans = 1; else base_pair[i].pair_to[k]->cis_or_trans = 2; } } for (l = 0; l < num; l++) { if (n_contact >= n_alloc) { n_alloc += 50; c_tmp = contact; contact = new CONTACT*[n_alloc]; for (m = 0; m < n_contact; m++) contact[m] = c_tmp[m]; for (m = n_contact; m < n_alloc; m++) contact[m] = NULL; delete [] c_tmp; } contact[n_contact] = new CONTACT; strcpy(contact[n_contact]->a_chnid, base[i]->res->GetFirstAtom()->chnid()); strcpy(contact[n_contact]->a_resnam, base[i]->res->ResName()); strcpy(contact[n_contact]->a_atomtyp, atom_pair[pos[l]].a_atomtyp); contact[n_contact]->a_atomnum = FstIndex; contact[n_contact]->a_resnum = FstRes; strcpy(contact[n_contact]->b_chnid, base[j]->res->GetFirstAtom()->chnid()); strcpy(contact[n_contact]->b_resnam, base[j]->res->ResName()); strcpy(contact[n_contact]->b_atomtyp, atom_pair[pos[l]].b_atomtyp); contact[n_contact]->b_atomnum = SndIndex; contact[n_contact]->b_resnum = SndRes; contact[n_contact]->type = base_pair[i].pair_to[k]->type; contact[n_contact]->sym = base[j]->op+1; contact[n_contact]->transl[0] = base[j]->lx; contact[n_contact]->transl[1] = base[j]->ly; contact[n_contact]->transl[2] = base[j]->lz; /* if (type_i && type_j && cis_or_trans) { sprintf(contact[n_contact]->eric_type, "%s %s/%s", _gb_dir[cis_or_trans], _edges[type_i], _edges[type_j]); } else strcpy(contact[n_contact]->eric_type, ""); */ n_contact++; } } } return contact; } static double cal_dist(_ATOM *atom_a, _ATOM *atom_b, const NDBCELLCONST *cell, const int op, const int lx, const int ly, const int lz) { double dx, dy, dz, dist, crysa, crysb, crysc; double cosalpha, cosbeta, cosgamma; if (!op && !lx && !ly && !lz) { dx = atom_a->orig().x - atom_b->orig().x; dy = atom_a->orig().y - atom_b->orig().y; dz = atom_a->orig().z - atom_b->orig().z; dist = sqrt(dx * dx + dy * dy + dz * dz); return dist; } COORD a = atom_a->orig(); COORD b = atom_b->orig(); ndb_trans_coord_3by3(a, cell, NDB_TRANS_ORTHO_TO_FRAC); ndb_trans_coord_3by3(b, cell, NDB_TRANS_ORTHO_TO_FRAC); do_symmetry_operation(b, b, &(cell->symops[op])); b.x += (double) lx; b.y += (double) ly; b.z += (double) lz; crysa = cell->a; crysb = cell->b; crysc = cell->c; cosalpha = cos(cell->alpha * degtorad); cosbeta = cos(cell->beta * degtorad); cosgamma = cos(cell->gamma * degtorad); dx = a.x - b.x; dy = a.y - b.y; dz = a.z - b.z; dist = sqrt(crysa * dx * crysa * dx + crysb * dy * crysb * dy + crysc * dz * crysc * dz + 2.0 * (crysa * dx * crysb * dy * cosgamma + crysa * dx * crysc * dz * cosbeta + crysb * dy * crysc * dz * cosalpha)); return dist; } static int get_best_pair(_BASE_PAIR *base_pair, int &pair_type, int &type_i, int &type_j, int &cis_or_trans, int *skip) { int i = 0, j = -1; double pair_score = XBIG; for (i = 0; i < base_pair->num; i++) { if (skip[base_pair->pair_to[i]->index]) continue; if (base_pair->pair_to[i]->type < 0) continue; if (base_pair->pair_to[i]->d2dp < pair_score) { /* if (base_pair->pair_to[i]->type == HBOND_TYPE_19 || base_pair->pair_to[i]->type == HBOND_TYPE_20) pair_type = 2; else pair_type = 0; */ pair_type = base_pair->pair_to[i]->type; type_i = base_pair->pair_to[i]->type_i; type_j = base_pair->pair_to[i]->type_j; cis_or_trans = base_pair->pair_to[i]->cis_or_trans; j = base_pair->pair_to[i]->index; pair_score = base_pair->pair_to[i]->d2dp; } } return j; } static void bpstep_par(COORD rot1[3], const COORD org1, COORD rot2[3], const COORD org2, double *pars, COORD mst_orien[3], COORD &mst_org) { double phi, rolltilt; COORD t1, hinge, mstx, msty, para_bp1[3], para_bp2[3], temp[3]; vector_cross_product(hinge, rot1[2], rot2[2]); rolltilt = magang(rot1[2], rot2[2]); arb_rotation(temp, hinge, -0.5 * rolltilt); multi_matrix(para_bp2, temp, rot2); arb_rotation(temp, hinge, 0.5 * rolltilt); multi_matrix(para_bp1, temp, rot1); pars[5] = vector_angle(para_bp1[1], para_bp2[1], para_bp2[2]); if (fabs(pars[5] - 180) < XEPS) { msty = para_bp2[0]; } else if (fabs(pars[6] + 180) < XEPS) { msty = para_bp2[0]; vector_reverse(msty); } else { vector_sum(msty, para_bp1[1], para_bp2[1]); vector_normalize(msty); } vector_cross_product(mstx, msty, para_bp2[2]); vector_middle(mst_org, org1, org2); vector_difference(t1, org2, org1); mst_orien[0] = mstx; mst_orien[1] = msty; mst_orien[2] = para_bp2[2]; pars[0] = vector_dot_product(t1, mstx); pars[1] = vector_dot_product(t1, msty); pars[2] = vector_dot_product(t1, para_bp2[2]); phi = deg2rad(vector_angle(hinge, msty, para_bp2[2])); pars[4] = rolltilt * cos(phi); pars[3] = rolltilt * sin(phi); } static void helical_par(COORD rot1[3], const COORD org1, COORD rot2[3], const COORD org2, double *pars, COORD mst_orien[3], COORD &mst_org) { double AD_mag, phi, TipInc1, TipInc2, vlen; COORD rot1_h[3], rot2_h[3], temp[3]; COORD t1, t2, axis_h, hinge1, hinge2, AD_axis, org1_h, org2_h; vector_difference(t1, rot2[0], rot1[0]); vector_difference(t2, rot2[1], rot1[1]); vector_cross_product(axis_h, t1, t2); vlen = vector_length(axis_h); if (vlen < XEPS) { axis_h.x = 0.0; axis_h.y = 0.0; axis_h.z = 1.0; } else { axis_h.x /= vlen; axis_h.y /= vlen; axis_h.z /= vlen; } TipInc1 = magang(axis_h, rot1[2]); vector_cross_product(hinge1, axis_h, rot1[2]); arb_rotation(temp, hinge1, -TipInc1); multi_matrix(rot1_h, temp, rot1); TipInc2 = magang(axis_h, rot2[2]); vector_cross_product(hinge2, axis_h, rot2[2]); arb_rotation(temp, hinge2, -TipInc2); multi_matrix(rot2_h, temp, rot2); vector_sum(t1, rot1_h[0], rot2_h[0]); vector_sum(t2, rot1_h[1], rot2_h[1]); vector_normalize(t1); vector_normalize(t2); mst_orien[0] = t1; mst_orien[1] = t2; mst_orien[2] = axis_h; pars[5] = vector_angle(rot1_h[1], rot2_h[1], axis_h); vector_difference(t2, org2, org1); pars[2] = vector_dot_product(t2, axis_h); phi = deg2rad(vector_angle(hinge1, rot1_h[1], axis_h)); pars[4] = TipInc1 * cos(phi); pars[3] = TipInc1 * sin(phi); vector_sum_scaled(t1, t2, -pars[2], axis_h); if (fabs(pars[5]) < HTWIST0) vector_sum_scaled(org1_h, org1, 0.5, t1); else { get_vector(t1, axis_h, 90 - pars[5] / 2.0, AD_axis); AD_mag = 0.5 * vector_length(t1) / sin(deg2rad(pars[5] / 2)); vector_sum_scaled(org1_h, org1, AD_mag, AD_axis); } vector_sum_scaled(org2_h, org1_h, pars[2], axis_h); vector_middle(mst_org, org1_h, org2_h); vector_difference(t1, org1, org1_h); pars[0] = vector_dot_product(t1, rot1_h[0]); pars[1] = vector_dot_product(t1, rot1_h[1]); } static void project_xyzP(const COORD P_1_i_plus_1, const COORD P_2_i, COORD mst_orien[3], COORD mst_org, COORD& aveP) { COORD temp, P_mst1, P_mst2; vector_difference(temp, P_1_i_plus_1, mst_org); P_mst1.x = vector_dot_product(temp, mst_orien[0]); P_mst1.y = vector_dot_product(temp, mst_orien[1]); P_mst1.z = vector_dot_product(temp, mst_orien[2]); vector_difference(temp, P_2_i, mst_org); P_mst2.x = vector_dot_product(temp, mst_orien[0]); P_mst2.y = -vector_dot_product(temp, mst_orien[1]); P_mst2.z = -vector_dot_product(temp, mst_orien[2]); vector_middle(aveP, P_mst1, P_mst2); } static double distance_ab(_BASE *a, const int ia, _BASE *b, const int ib) { COORD pa, pb, pc; if (ia) pa = a->P_xyz; else pa = a->O_xyz; if (ib) pb = b->P_xyz; else pb = b->O_xyz; vector_difference(pc, pa, pb); return (vector_length(pc)); } static void get_ij(const int m, const int index, int *best_pairs, int &i, int &j) { i = best_pairs[6 * m]; j = best_pairs[6*m+1]; if (index) { i = best_pairs[6*m+1]; j = best_pairs[6 * m]; } } static void reverse(const int st, const int n, int *vec) { int *vtmp, i = 0; vtmp = new int[n]; for (i = 0; i < n; i++) vtmp[i] = vec[n - 1 + st - i]; for (i = 0; i < n; i++) vec[st + i] = vtmp[i]; delete [] vtmp; } static int is_linked(_BASE *base1, _BASE *base2) { if (base1->exist_O && base2->exist_P) { if (distance_ab(base1, 0, base2, 1) <= O3P_UPPER) return 1; } if (base1->exist_P && base2->exist_O) { if (distance_ab(base1, 1, base2, 0) <= O3P_UPPER) return -1; } return 0; } static int is_wc_geometry(_BASE *base1, _BASE *base2) { double dir_x = vector_dot_product(base1->orien[0], base2->orien[0]); double dir_y = vector_dot_product(base1->orien[1], base2->orien[1]); double dir_z = vector_dot_product(base1->orien[2], base2->orien[2]); return (dir_x > 0.0 && dir_y < 0.0 && dir_z < 0.0); } static int check_others(_BASE **_base, const int i1, const int j1, const int i2, const int j2) { int ii1, ii2, jj1, jj2; double a1[3], a2[3], r1[3], r2[3], d[4]; if (is_linked(_base[i1], _base[i2]) || is_linked(_base[j1], _base[j2]) || is_linked(_base[i1], _base[j2]) || is_linked(_base[j1], _base[i2])) return 0; a1[0] = vector_dot_product(_base[i1]->orien[0], _base[i2]->orien[0]); a1[1] = vector_dot_product(_base[i1]->orien[1], _base[i2]->orien[1]); a1[2] = vector_dot_product(_base[i1]->orien[2], _base[i2]->orien[2]); a2[0] = vector_dot_product(_base[j1]->orien[0], _base[j2]->orien[0]); a2[1] = vector_dot_product(_base[j1]->orien[1], _base[j2]->orien[1]); a2[2] = vector_dot_product(_base[j1]->orien[2], _base[j2]->orien[2]); ii1 = a1[0] > 0 && a1[1] > 0 && a1[2] > 0; ii2 = a2[0] > 0 && a2[1] > 0 && a2[2] > 0; if (ii1 && ii2) return 0; r1[0] = vector_dot_product(_base[i1]->orien[0], _base[j2]->orien[0]); r1[1] = vector_dot_product(_base[i1]->orien[1], _base[j2]->orien[1]); r1[2] = vector_dot_product(_base[i1]->orien[2], _base[j2]->orien[2]); r2[0] = vector_dot_product(_base[j1]->orien[0], _base[i2]->orien[0]); r2[1] = vector_dot_product(_base[j1]->orien[1], _base[i2]->orien[1]); r2[2] = vector_dot_product(_base[j1]->orien[2], _base[i2]->orien[2]); jj1 = r1[0] > 0 && r1[1] > 0 && r1[2] > 0; jj2 = r2[0] > 0 && r2[1] > 0 && r2[2] > 0; if (!ii1 && !ii2) { if (jj1 || jj2) return 1; else if (!jj1 && !jj2) return 0; } d[0] = dot2ang(a1[0]) + dot2ang(a1[1]) + dot2ang(a1[2]); d[1] = dot2ang(a2[0]) + dot2ang(a2[1]) + dot2ang(a2[2]); d[2] = dot2ang(r1[0]) + dot2ang(r1[1]) + dot2ang(r1[2]); d[3] = dot2ang(r2[0]) + dot2ang(r2[1]) + dot2ang(r2[2]); if (ii1 && jj1) return (d[0] > d[2]) ? 1 : 0; if (ii1 && jj2) return (d[0] > d[3]) ? 1 : 0; if (ii2 && jj1) return (d[1] > d[2]) ? 1 : 0; if (ii2 && jj2) return (d[1] > d[3]) ? 1 : 0; return 0; } static int check_direction(int *helix_idx, int *bp_idx, int *index, int *best_pairs, _BASE **_base, int *direction) { int i1, j1, i2, j2, j, k, m, n; for (j = 0; j < 6; j++) direction[j] = 0; for (j = helix_idx[0]; j < helix_idx[1] - 1; j++) { m = bp_idx[j]; n = bp_idx[j+1]; get_ij(m, index[m], best_pairs, i1, j1); get_ij(n, index[n], best_pairs, i2, j2); k = is_linked(_base[i1], _base[i2]); (k == 1) ? ++direction[0] : (k == -1) ? ++direction[1] : ++direction[2]; k = is_linked(_base[j1], _base[j2]); (k == 1) ? ++direction[3] : (k == -1) ? ++direction[4] : ++direction[5]; } if ((direction[0] && direction[1]) || (direction[3] && direction[4])) return 0; if ((direction[0] + direction[1] + direction[3] + direction[4]) == 0) return 1; m = bp_idx[helix_idx[0]]; n = bp_idx[helix_idx[1] - 1]; get_ij(m, index[m], best_pairs, i1, j1); get_ij(n, index[n], best_pairs, i2, j2); if (direction[0] && !direction[1]) { if (!direction[3] && direction[4]) { if (i1 > j2) { for (j = helix_idx[0]; j < helix_idx[1]; j++) index[bp_idx[j]] = !index[bp_idx[j]]; reverse(helix_idx[0], helix_idx[1] - helix_idx[0], bp_idx); } } else if (direction[3] && !direction[4]) { if (i1 > j1) { for (j = helix_idx[0]; j < helix_idx[1]; j++) index[bp_idx[j]] = !index[bp_idx[j]]; } } } return 1; } static int find_atom_index(const int num, _ATOMNAME *name_list, const char *name) { int index = -1; _ATOMNAME tmp, *exist = NULL; strcpy(tmp.name, name); exist = (_ATOMNAME *) bsearch(&tmp, name_list, num, sizeof(_ATOMNAME), name_compare1); if (exist) index = ((long int) exist - (long int) name_list) / sizeof(_ATOMNAME); return index; } #ifdef SUN_OS #ifdef __cplusplus extern "C" { #endif #endif static int base_hbond_type_compare(const void *base_type1, const void *base_type2) { return (strcmp(((BASE_HBOND_TYPE *) base_type1)->type, ((BASE_HBOND_TYPE *) base_type2)->type)); } static int name_compare1(const void *name1, const void *name2) { return (strcmp(((_ATOMNAME *) name1)->name, ((_ATOMNAME *) name2)->name)); } static int base_type_compare(const void *base_type1, const void *base_type2) { return (strcmp(((_DONOR_ACCEPTOR *) base_type1)->type, ((_DONOR_ACCEPTOR *) base_type2)->type)); } static int atom_edge_compare(const void *atom_edge1, const void *atom_edge2) { return (strcmp(((_Atom_Edge *) atom_edge1)->atom, ((_Atom_Edge *) atom_edge2)->atom)); } #ifdef SUN_OS #ifdef __cplusplus } #endif #endif