forked from pemsley/coot
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsimple-restraint.cc
5030 lines (4190 loc) · 172 KB
/
simple-restraint.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* ideal/simple-restraint.cc
*
* Copyright 2002, 2003, 2004, 2005, 2006 by The University of York
* Copyright 2008, 2009, 2010 by The University of Oxford
* Copyright 2013, 2014, 2015, 2016 by Medical Research Council
* Author: Paul Emsley
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU General Public License and
* the GNU Lesser General Public License along with this program; if not,
* write to the Free Software Foundation, Inc., 51 Franklin Street,
* Fifth Floor, Boston, MA, 02110-1301, USA.
*/
// #define ANALYSE_REFINEMENT_TIMING
#include <string.h> // for strcmp
#ifdef ANALYSE_REFINEMENT_TIMING
#ifdef _MSC_VER
#include <time.h>
#else
#include <sys/time.h>
#endif
#endif // ANALYSE_REFINEMENT_TIMING
#include <fstream>
#include <algorithm> // for sort
#include <stdexcept>
#include <iomanip>
#ifdef HAVE_CXX_THREAD
#include <thread>
#include <chrono>
#endif // HAVE_CXX_THREAD
#include "utils/split-indices.hh"
#include "geometry/mol-utils.hh"
#include "geometry/main-chain.hh"
#include "simple-restraint.hh"
//
#include "coot-utils/coot-coord-extras.hh" // is_nucleotide_by_dict
// #include "mmdb.h" // for printing of mmdb::Atom pointers as info not raw
// pointers. Removed. Too much (linking issues in)
// Makefile pain.
#include "compat/coot-sysdep.h"
zo::rama_table_set coot::restraints_container_t::zo_rama;
std::atomic<bool> coot::restraints_container_t::print_lock(false);
void
coot::restraints_container_t::clear() {
bool unlocked = false;
while (! restraints_lock.compare_exchange_weak(unlocked, true)) {
std::this_thread::sleep_for(std::chrono::microseconds(10));
unlocked = false;
}
restraints_vec.clear();
init(); // resets lock, fwiw
}
void
coot::restraints_container_t::get_restraints_lock() {
bool unlocked = false;
while (! restraints_lock.compare_exchange_weak(unlocked, true)) {
std::this_thread::sleep_for(std::chrono::nanoseconds(10));
unlocked = false;
}
}
void
coot::restraints_container_t::release_restraints_lock() {
restraints_lock = false;
}
// static
void
coot::restraints_container_t::get_print_lock() {
bool unlocked = false;
while (! print_lock.compare_exchange_weak(unlocked, true)) {
std::this_thread::sleep_for(std::chrono::microseconds(1));
unlocked = false;
}
}
// static
void
coot::restraints_container_t::release_print_lock() {
print_lock = false;
}
coot::restraints_container_t::~restraints_container_t() {
unset_fixed_during_refinement_udd();
// 20240228-PE If the restraints_container_t is closed/finished with before the refinement
// has terminated (i.e. continue_status == GSL_CONTINUE), then free_delete_reset()
// doesn't get called. So call it now.
// Fixes memory leak.
//
// 20240229-PE Nope. This causes a crash is test_peptide_omega()
// free_delete_reset();
if (from_residue_vector) {
if (atom) {
// this is constructed manually.
// Oh we can't do this here because we copy the
// restraints in simple_refine_residues() and that
// shallow copies the atom pointer - the original
// restraints go out of scope and call this destructor.
//
// We need a new way to get rid of atom - c.f. the
// linear/conventional way?
//
// 20200820-PE I don't know what simple_refine_residues() is
// I am going to ignore the above message and delete the atoms
// now.
// atom needs to be a shared_ptr so that I can copy
// restraints containers.
// delete [] atom;
// atom = NULL;
// 20240228-PE another memory leak: atom
// It seems to me at the moment that atom can be assigned in the constructor
// in which case, we don't want to delete it.
// Or it can be assigned in init_from_residue_vec(), in which case, we do
// want to delete it (atom_array_needs_to_be_deleted_at_end is set there)
if (atom_array_needs_to_be_deleted_at_end) {
delete [] atom;
atom = nullptr;
}
}
} else {
// member data item mmdb::PPAtom atom is constructed by an
// mmdb SelectAtoms()/GetSelIndex() (which includes
// flanking atoms).
// 20081207: don't do this here now - because the
// memory/selection is deleted again in
// clear_up_moving_atoms(). It *should* be done here of
// course, but we'll save that for the future.
//
//if (atom) {
// mol->DeleteSelection(SelHnd_atom);
// atom = NULL;
// }
}
}
// Used in omega distortion graph
//
coot::restraints_container_t::restraints_container_t(atom_selection_container_t asc_in,
const std::string &chain_id,
const clipper::Xmap<float> *map_p_in) : xmap_p(map_p_in) {
init();
mol = asc_in.mol;
are_all_one_atom_residues = false;
istart_res = 999999;
iend_res = -9999999;
mmdb::PResidue *SelResidues = NULL;
int nSelResidues;
// -------- Find the max and min res no -----------------------------
int selHnd = mol->NewSelection();
mol->Select(selHnd, mmdb::STYPE_RESIDUE, 1,
chain_id.c_str(),
mmdb::ANY_RES, "*",
mmdb::ANY_RES, "*",
"*", // residue name
"*", // Residue must contain this atom name?
"*", // Residue must contain this Element?
"*", // altLocs
mmdb::SKEY_NEW // selection key
);
mol->GetSelIndex(selHnd, SelResidues, nSelResidues);
for (int ires=0; ires<nSelResidues; ires++) {
int resno = SelResidues[ires]->GetSeqNum();
if (resno < istart_res)
istart_res = resno;
if (resno > iend_res)
iend_res = resno;
}
mol->DeleteSelection(selHnd);
//
// -------- Found the max and min res no -----------------------------
// -------------------------------------------------------------------
// Set class variables atom to the selection that includes the
// chain (which is not the same as the input atom selection)
//
int SelHnd = mol->NewSelection();
atom = NULL;
mol->SelectAtoms(SelHnd, 0,
chain_id.c_str(),
mmdb::ANY_RES, // starting resno, an int
"*", // any insertion code
mmdb::ANY_RES, // ending resno
"*", // ending insertion code
"*", // any residue name
"*", // atom name
"*", // elements
"*" // alt loc.
);
mol->GetSelIndex(SelHnd, atom, n_atoms);
// -------------------------------------------------------------------
initial_position_params_vec.resize(3*n_atoms);
for (int i=0; i<n_atoms; i++) {
initial_position_params_vec[3*i ] = atom[i]->x;
initial_position_params_vec[3*i+1] = atom[i]->y;
initial_position_params_vec[3*i+2] = atom[i]->z;
// std::cout << " " << i << " " << coot::atom_spec_t(atom[i]) << "\n";
}
}
coot::restraints_container_t::restraints_container_t(mmdb::PResidue *SelResidues, int nSelResidues,
const std::string &chain_id,
mmdb::Manager *mol_in,
const clipper::Xmap<float> *map_p_in) : xmap_p(map_p_in) {
init();
are_all_one_atom_residues = false;
std::vector<coot::atom_spec_t> fixed_atoms_dummy;
int istart_res_l = 999999;
int iend_res_l = -9999999;
int resno;
for (int i=0; i<nSelResidues; i++) {
resno = SelResidues[i]->seqNum;
if (resno < istart_res_l)
istart_res_l = resno;
if (resno > iend_res_l)
iend_res_l = resno;
}
short int have_flanking_residue_at_start = 0;
short int have_flanking_residue_at_end = 0;
short int have_disulfide_residues = 0;
const char *chn = chain_id.c_str();
// std::cout << "DEBUG: ==== istart_res iend_res " << istart_res << " "
// << iend_res << std::endl;
init_from_mol(istart_res_l, iend_res_l,
have_flanking_residue_at_start,
have_flanking_residue_at_end,
have_disulfide_residues,
std::string(""), chn, mol_in, fixed_atoms_dummy);
}
bool
coot::residue_sorter(const std::pair<bool, mmdb::Residue *> &r1,
const std::pair<bool, mmdb::Residue *> &r2) {
std::string chain_id_1 = r1.second->GetChainID();
std::string chain_id_2 = r2.second->GetChainID();
if (chain_id_1 < chain_id_2) {
return true;
} else {
if (chain_id_1 > chain_id_2) {
return false;
} else {
if (r1.second->index < r2.second->index) {
return true;
} else {
if (r1.second->index > r2.second->index) {
return false;
} else {
if (r1.second->GetSeqNum() < r2.second->GetSeqNum()) {
return true;
} else {
if (r1.second->GetSeqNum() > r2.second->GetSeqNum()) {
return false;
} else {
std::string ins_code_1 = r1.second->GetInsCode();
std::string ins_code_2 = r2.second->GetInsCode();
if (ins_code_1 < ins_code_2) {
return true;
} else {
if (ins_code_1 > ins_code_2) {
return false;
}
}
}
}
}
}
}
}
return false;
}
// 20081106 construct from a vector of residues, each of which
// has a flag attached that denotes whether or not it is a fixed
// residue (it would be set, for example in the case of flanking
// residues).
//
// currently links are ignored.
coot::restraints_container_t::restraints_container_t(const std::vector<std::pair<bool,mmdb::Residue *> > &residues,
const std::vector<mmdb::Link> &links_in,
const coot::protein_geometry &geom,
mmdb::Manager *mol_in,
const std::vector<atom_spec_t> &fixed_atom_specs,
const clipper::Xmap<float> *map_p_in) : xmap_p(map_p_in) {
istart_minus_flag = false; // used in make_flanking_atoms_rama_restraints
iend_plus_flag = false;
init();
from_residue_vector = 1;
are_all_one_atom_residues = false;
// int n_resiudes_in_mol = util::number_of_residues_in_molecule(mol_in);
// std::cout << "################## n_resiudes_in_mol " << n_resiudes_in_mol << std::endl;
std::vector<std::pair<bool,mmdb::Residue *> > residues_local;
residues_local.reserve(residues.size());
for(unsigned int i=0; i<residues.size(); i++)
if (residues[i].second)
residues_local.push_back(residues[i]);
// now sort those residues so that polymer linking is easy
// sorting function should return false at the end - because sorting function is called
// with the same argument for left and right hand side, I think (std::sort testing for sanity?)
std::sort(residues_local.begin(), residues_local.end(), residue_sorter);
if (false)
for (std::size_t i=0; i<residues_local.size(); i++)
std::cout << "DEBUG:: restraints_container_t() constructor: passed-residue-vec residue: "
<< residue_spec_t(residues_local[i].second)
<< " has index " << residues_local[i].second->index << std::endl;
residues_vec = residues_local;
init_from_residue_vec(residues_local, geom, mol_in, fixed_atom_specs);
fill_links(mol_in);
}
coot::restraints_container_t::restraints_container_t(const std::vector<std::pair<bool,mmdb::Residue *> > &residues,
const coot::protein_geometry &geom,
mmdb::Manager *mol_in,
const clipper::Xmap<float> *map_p_in) : xmap_p(map_p_in) {
init();
from_residue_vector = 1;
are_all_one_atom_residues = false;
std::vector<atom_spec_t> dummy_fixed_atom_specs;
std::vector<std::pair<bool,mmdb::Residue *> > residues_local;
residues_local.reserve(residues.size());
for(unsigned int i=0; i<residues.size(); i++)
if (residues[i].second)
residues_local.push_back(residues[i]);
std::sort(residues_local.begin(), residues_local.end(), residue_sorter);
residues_vec = residues_local;
if (false) {
std::cout << "debug:: in restraints_container_t() constructor with local residue size " << residues_local.size() << std::endl;
for (unsigned int ir=0; ir<residues_vec.size(); ir++) {
if (residues_vec[ir].second) {
std::cout << "INFO:: starting init_from_residue_vec() residue " << residues_vec[ir].second << std::endl;
} else {
std::cout << "ERROR:: starting init_from_residue_vec() NUll residue " << ir << std::endl;
}
}
}
//std::cout << "##################################### constructor D "
init_from_residue_vec(residues_local, geom, mol_in, dummy_fixed_atom_specs); // sets mol
fill_links(mol_in);
}
// What are the rules for dealing with alt conf in flanking residues?
//
// First we want to try to select only atom that have the same alt
// conf, if that fails to select atoms in flanking residues (and there
// are atoms in flanking residues (which we know from
// have_flanking_residue_at_end/start), then we should try an mmdb construction [""|"A"]
// (i.e. blank or "A"). If that fails, give up - it's a badly formed pdb file (it
// seems to me).
//
void
coot::restraints_container_t::init_from_mol(int istart_res_in, int iend_res_in,
bool have_flanking_residue_at_start,
bool have_flanking_residue_at_end,
short int have_disulfide_residues,
const std::string &altloc,
const std::string &chain_id,
mmdb::Manager *mol_in,
const std::vector<coot::atom_spec_t> &fixed_atom_specs) {
init_shared_pre(mol_in);
istart_res = istart_res_in;
iend_res = iend_res_in;
chain_id_save = chain_id;
// internal flags that mirror passed have_flanking_residue_at_* variables
istart_minus_flag = have_flanking_residue_at_start;
iend_plus_flag = have_flanking_residue_at_end;
int iselection_start_res = istart_res;
int iselection_end_res = iend_res;
// std::cout << "start res range: " << istart_res << " " << iend_res << " " << chain_id << "\n";
// Are the flanking atoms available in mol_in?
// mol_in was constructed outside of this class, so the mol_in constructing
// routine knows if they were there or not.
//
if (have_flanking_residue_at_start) iselection_start_res--;
if (have_flanking_residue_at_end) iselection_end_res++;
SelHnd_atom = mol->NewSelection();
mol->SelectAtoms(SelHnd_atom,
0,
chain_id.c_str(),
iselection_start_res, "*",
iselection_end_res, "*",
"*", // rnames
"*", // anames
"*", // elements
"*" // altLocs
);
// set the mmdb::PPAtom atom (class variable) and n_atoms:
//
mol->GetSelIndex(SelHnd_atom, atom, n_atoms);
if (false) { // debugging;
std::cout << "debug:: in init_from_mol() here are the " << fixed_atom_indices.size()
<< " fixed_atom indices: \n";
std::set<int>::const_iterator it;
for (it=fixed_atom_indices.begin(); it!=fixed_atom_indices.end(); ++it)
std::cout << " " << *it;
std::cout << "\n";
for (int iat=0; iat<n_atoms; iat++)
std::cout << " " << iat << " " << coot::atom_spec_t(atom[iat]) << " with altloc :"
<< altloc << ":" << std::endl;
}
bool debug = false;
if (debug) {
std::cout << "DEBUG:: Selecting residues in chain \"" << chain_id << "\" gives "
<< n_atoms << " atoms " << std::endl;
for (int iat=0; iat<n_atoms; iat++) {
std::cout << " " << iat << " " << atom[iat]->name << " " << atom[iat]->GetSeqNum()
<< " " << atom[iat]->GetChainID() << std::endl;
}
}
// debugging, you need to uncomment mmdb.h at the top and add coords to the linking
// atom_selection_container_t tmp_res_asc = make_asc(mol_in);
// std::cout << "There are " << tmp_res_asc.n_selected_atoms
// << " atoms in tmp_res_asc\n";
// for (int kk=0; kk<tmp_res_asc.n_selected_atoms; kk++) {
// std::cout << "In simple rest " << kk << " "
// << tmp_res_asc.atom_selection[kk] << "\n";
// }
// std::cout << "INFO::" << n_atoms << " atoms selected from molecule for refinement" ;
// std::cout << " (this includes fixed and flanking atoms)." << std::endl;
if (n_atoms == 0) {
std::cout << "ERROR:: atom selection disaster:" << std::endl;
std::cout << " This should not happen" << std::endl;
std::cout << " residue range: " << iselection_start_res << " "
<< iselection_end_res << " chain-id \"" << chain_id << "\" "
<< "flanking flags: " << have_flanking_residue_at_start
<< " " << have_flanking_residue_at_end << std::endl;
}
init_shared_post(fixed_atom_specs); // clears fixed_atom_indices
add_fixed_atoms_from_flanking_residues(have_flanking_residue_at_start,
have_flanking_residue_at_end,
iselection_start_res, iselection_end_res);
}
void
coot::restraints_container_t::init_shared_pre(mmdb::Manager *mol_in) {
needs_reset = false;
verbose_geometry_reporting = NORMAL;
do_numerical_gradients_flag = false;
have_oxt_flag = false; // set in mark_OXT()
dist_crit_for_bonded_pairs = 3.0;
// the smaller the alpha, the more like least squares
geman_mcclure_alpha = 0.2; // Is this a good value? Talk to Rob.
mol = mol_in;
lennard_jones_epsilon = 1.0; // close enough to 0.997 kJ/mol
cryo_em_mode = true;
n_times_called = 0;
n_small_cycles_accumulator = 0;
m_s = 0;
x = 0;
#ifdef HAVE_BOOST_BASED_THREAD_POOL_LIBRARY
n_threads = 0;
#endif // HAVE_BOOST_BASED_THREAD_POOL_LIBRARY
log_cosh_target_distance_scale_factor = 3000.0;
convert_plane_restraints_to_improper_dihedral_restraints_flag = false; // as it was in 2019.
use_proportional_editing = false;
pull_restraint_neighbour_displacement_max_radius = 10.0; // make this a member of the class
init_neutron_occupancies();
}
double
coot::restraints_container_t::get_distortion_score() const {
// Yummy mixing C and C++ APIs...
return distortion_score(x, const_cast<void *>(static_cast<const void *>(this)));
}
void
coot::restraints_container_t::set_use_proportional_editing(bool state) {
use_proportional_editing = state;
}
void
coot::restraints_container_t::set_has_hydrogen_atoms_state() {
// in init model_has_hydrogens is set to true;
bool found = false;
for (int i=0; i<n_atoms; i++) {
mmdb::Atom *at = atom[i];
if (is_hydrogen(at)) {
found = true;
break;
}
}
if (! found)
model_has_hydrogen_atoms = false;
}
// pass the formal charge also?
double
coot::restraints_container_t::neutron_occupancy(const std::string &element, int formal_charge) const {
std::string mod_ele = coot::util::remove_whitespace(element);
if (mod_ele.length() > 1)
mod_ele = coot::util::capitalise(mod_ele);
if (formal_charge != 0)
mod_ele += coot::util::int_to_string(formal_charge);
std::map<std::string, double>::const_iterator it = neutron_occupancy_map.find(mod_ele);
if (it != neutron_occupancy_map.end())
return it->second;
else
return 0.0;
}
void
coot::restraints_container_t::set_z_occ_weights() {
// z weights:
//
atom_z_occ_weight.resize(n_atoms);
std::vector<std::pair<std::string, int> > atom_list = coot::util::atomic_number_atom_list();
for (int i=0; i<n_atoms; i++) {
mmdb::Atom *at = atom[i];
if (! at->isTer()) {
std::string element = at->element;
double z = coot::util::atomic_number(at->element, atom_list);
double weight = 1.0;
double occupancy = atom[i]->occupancy;
if (occupancy > 1.0) occupancy = 1.0;
if (do_neutron_refinement) {
int formal_charge = 0;
occupancy = neutron_occupancy(element, formal_charge);
}
if (cryo_em_mode) {
// is-side-chain? would be a better test
if (! is_main_chain_or_cb_p(at)) {
// std::cout << "downweighting atom " << coot::atom_spec_t(atom[i]) << std::endl;
weight = 0.2;
}
std::string at_name = atom[i]->name;
if (at_name == " O ") {
weight = 0.4;
}
}
if (z < 0.0) {
std::cout << "WARNING:: init_shared_post() atom " << i << " " << atom_spec_t(atom[i])
<< " Unknown element \"" << atom[i]->element << "\"" << std::endl;
z = 6.0; // as for carbon
}
atom_z_occ_weight[i] = weight * z * occupancy;
}
}
}
void
coot::restraints_container_t::init_shared_post(const std::vector<atom_spec_t> &fixed_atom_specs) {
// std::cout << "##################################### init_shared_post() " << fixed_atom_specs.size() << std::endl;
bonded_atom_indices.resize(n_atoms);
set_has_hydrogen_atoms_state();
initial_position_params_vec.resize(3*n_atoms);
for (int i=0; i<n_atoms; i++) {
initial_position_params_vec[3*i ] = atom[i]->x;
initial_position_params_vec[3*i+1] = atom[i]->y;
initial_position_params_vec[3*i+2] = atom[i]->z;
}
// Set the UDD have_bond_or_angle to initally all "not". They get
// set to "have" (1) in make_restraints (and functions thereof).
//
// udd_handle becomes member data so that it can be used in
// make_restraints() without passing it back (this function is part
// of a constructor, don't forget).
//
// 20131213-PE: I dont see the point of udd_bond_angle.
//
if (mol) {
udd_bond_angle = mol->RegisterUDInteger (mmdb::UDR_ATOM, "bond or angle");
if (udd_bond_angle < 0) {
std::cout << "ERROR:: can't make udd_handle in init_from_mol\n";
} else {
for (int i=0; i<n_atoms; i++) {
atom[i]->PutUDData(udd_bond_angle,0);
}
}
}
// Set the UDD of the indices in the atom array (i.e. the thing
// that get_asc_index returns)
//
if (mol) {
udd_atom_index_handle = mol->RegisterUDInteger ( mmdb::UDR_ATOM, "atom_array_index");
if (udd_atom_index_handle < 0) {
std::cout << "ERROR:: can't make udd_handle in init_from_mol\n";
} else {
for (int i=0; i<n_atoms; i++) {
atom[i]->PutUDData(udd_atom_index_handle,i);
// std::cout << "init_shared_post() atom " << atom_spec_t(atom[i])
// << " gets udd_atom_index_handle value " << i << std::endl;
}
}
}
use_map_gradient_for_atom.resize(n_atoms, false);
if (! from_residue_vector) {
// convential way
for (int i=0; i<n_atoms; i++) {
if (atom[i]->residue->seqNum >= istart_res &&
atom[i]->residue->seqNum <= iend_res) {
if (! is_hydrogen(atom[i]))
use_map_gradient_for_atom[i] = true;
} else {
use_map_gradient_for_atom[i] = false;
}
}
} else {
// blank out the non moving atoms (i.e. flanking residues)
for (int i=0; i<n_atoms; i++) {
mmdb::Residue *res_p = atom[i]->residue;
if (is_a_moving_residue_p(res_p)) {
if (! is_hydrogen(atom[i]) || do_hydrogen_atom_refinement)
use_map_gradient_for_atom[i] = true;
} else {
// std::cout << "blanking out density for atom " << i << std::endl;
use_map_gradient_for_atom[i] = false;
}
}
}
set_z_occ_weights();
// the fixed atoms:
//
assign_fixed_atom_indices(fixed_atom_specs); // convert from std::vector<atom_spec_t>
// to std::vector<int> fixed_atom_indices;
// blank out those atoms from seeing electron density map gradients
std::set<int>::const_iterator it;
for (it=fixed_atom_indices.begin(); it!=fixed_atom_indices.end(); ++it)
use_map_gradient_for_atom[*it] = false;
if (verbose_geometry_reporting == VERBOSE)
for (int i=0; i<n_atoms; i++)
std::cout << atom[i]->name << " " << atom[i]->residue->seqNum << " "
<< use_map_gradient_for_atom[i] << std::endl;
}
// uses fixed_atom_indices
void
coot::restraints_container_t::set_fixed_during_refinement_udd() {
if (! mol) {
std::cout << "ERROR:: in set_fixed_during_refinement_udd() mol is null" << std::endl;
return;
}
int uddHnd = mol->RegisterUDInteger(mmdb::UDR_ATOM , "FixedDuringRefinement");
for (int i=0; i<n_atoms; i++) {
mmdb::Atom *at = atom[i];
//std::cout << " setting fixed udd flag on atom " << atom_spec_t(at) << " residue "
// << at->residue << " atom " << at << " mol " << mol << std::endl;
// if (std::find(fixed_atom_indices.begin(), fixed_atom_indices.end(), i) == fixed_atom_indices.end())
if (fixed_atom_indices.find(i) == fixed_atom_indices.end())
at->PutUDData(uddHnd, 0);
else
at->PutUDData(uddHnd, 1);
}
}
void
coot::restraints_container_t::unset_fixed_during_refinement_udd() {
if (! mol) {
std::cout << "ERROR:: in unset_fixed_during_refinement_udd() mol is null" << std::endl;
return;
}
int uddHnd = mol->GetUDDHandle(mmdb::UDR_ATOM , "FixedDuringRefinement");
for (int i=0; i<n_atoms; i++) {
mmdb::Atom *at = atom[i];
at->PutUDData(uddHnd, 0);
}
}
void
coot::restraints_container_t::init_from_residue_vec(const std::vector<std::pair<bool,mmdb::Residue *> > &residues,
const coot::protein_geometry &geom,
mmdb::Manager *mol_in,
const std::vector<atom_spec_t> &fixed_atom_specs) {
// std::cout << "##################################### start of init_from_residue_vec() "
// << fixed_atom_specs.size() << std::endl;
// This function is called from the constructor.
// make_restraints() is called after this function by the user of this class.
init_shared_pre(mol_in);
if (residues.size() > 2000000) {
std::cout << "ERROR:: in init_from_residue_vec() - memory error " << residues.size() << std::endl;
return;
}
residues_vec.resize(residues.size());
if (false) {
for (unsigned int ir=0; ir<residues_vec.size(); ir++) {
if (residues_vec[ir].second) {
std::cout << "INFO:: starting init_from_residue_vec() residue " << residues_vec[ir].second << std::endl;
} else {
std::cout << "ERROR:: starting init_from_residue_vec() NUll residue " << ir << std::endl;
}
}
}
// residues_vec = residues;
for (std::size_t i=0; i<residues.size(); i++)
if (residues[i].first == false) // i.e. is moving
residues_vec_moving_set.insert(residues[i].second);
// Need to set class members mmdb::PPAtom atom and int n_atoms.
// ...
// 20090620: or do we?
// debug:
bool debug = false;
if (debug) {
for (unsigned int ir=0; ir<residues_vec.size(); ir++) {
mmdb::PAtom *res_atom_selection = NULL;
int n_res_atoms;
residues_vec[ir].second->GetAtomTable(res_atom_selection, n_res_atoms);
std::cout << "debug:: =============== in init_from_residue_vec() residue "
<< ir << " of " << residues_vec.size() << " has : "
<< n_res_atoms << " atom " << std::endl;
std::cout << "debug:: =============== in init_from_residue_vec() residue "
<< ir << " of " << residues_vec.size() << " " << residue_spec_t(residues_vec[ir].second)
<< std::endl;
if (false)
for (int iat=0; iat<n_res_atoms; iat++) {
mmdb::Atom *at = res_atom_selection[iat];
std::cout << "DEBUG:: in init_from_residue_vec: atom "
<< iat << " of " << n_res_atoms << " \""
<< at->name << "\" \"" << at->altLoc << "\" "
<< at->GetSeqNum() << " \"" << at->GetInsCode() << "\" \""
<< at->GetChainID() << "\"" << std::endl;
}
}
}
// what about adding the flanking residues? How does the atom
// indexing of that work when (say) adding a bond?
if (false)
std::cout << "debug::info in init_from_residue_vec() calling bonded_flanking_residues_by_residue_vector() "
<< std::endl;
float dist_crit = 5.3; // 20170924-PE was 3.0 but this made a horrible link in a tight turn
// (which I suspect is not uncommon) crazy-neighbour-refine-519.pdb
// for EMDB 6224.
// 520 was bonded to 522 in a neighb (3-residue) refine on 519.
// This function is called by init but not make_restraints.
// init doesn't set bonded_pairs_container (make_restraints does that).
// dist_crit = 8.0
// fill neighbour_set from rnr (excluding residues of residues_vec):
std::vector<std::pair<bool,mmdb::Residue *> > active_residues_vec;
active_residues_vec.reserve(residues_vec.size()/2 + 1);
for (unsigned int ir=0; ir<residues_vec.size(); ir++) {
if (!residues_vec[ir].first)
active_residues_vec.push_back(residues_vec[ir]);
}
std::map<mmdb::Residue *, std::set<mmdb::Residue *> > rnr = residues_near_residues(active_residues_vec, mol, dist_crit);
std::map<mmdb::Residue *, std::set<mmdb::Residue *> > neighbour_set;
std::map<mmdb::Residue *, std::set<mmdb::Residue *> >::const_iterator it_map;
for(it_map=rnr.begin(); it_map!=rnr.end(); ++it_map) {
mmdb::Residue *r = it_map->first;
const std::set<mmdb::Residue *> &s = it_map->second;
std::set<mmdb::Residue *>::const_iterator it_set;
for (it_set=s.begin(); it_set!=s.end(); ++it_set) {
bool found = false;
for (std::size_t i=0; i<residues_vec.size(); i++) {
if (*it_set == residues_vec[i].second) {
if (! residues_vec[i].first) { // moving residue
found = true;
break;
}
}
}
if (! found) {
neighbour_set[r].insert(*it_set);
}
}
}
bonded_pair_container_t bpc = bonded_flanking_residues_by_residue_vector(neighbour_set, geom);
// std::cout << "bonded_pair_container_t bpc for flanking: " << std::endl;
// std::cout << bpc << std::endl;
// internal variable non_bonded_neighbour_residues is set by this
// function:
set_non_bonded_neighbour_residues_by_residue_vector(neighbour_set, bpc, geom);
if (false) { // debug
std::cout << "############## rnr map set: " << std::endl;
for(it_map=rnr.begin(); it_map!=rnr.end(); ++it_map) {
const std::set<mmdb::Residue *> &s = it_map->second;
std::cout << "::: rnr " << residue_spec_t(it_map->first) << std::endl;
std::set<mmdb::Residue *>::const_iterator it_set;
for (it_set=s.begin(); it_set!=s.end(); ++it_set) {
std::cout << " rnr " << residue_spec_t(*it_set) << std::endl;
}
}
std::cout << "############## neighbour set: " << std::endl;
for(it_map=neighbour_set.begin(); it_map!=neighbour_set.end(); ++it_map) {
const std::set<mmdb::Residue *> &s = it_map->second;
std::cout << "::: " << residue_spec_t(it_map->first) << std::endl;
std::set<mmdb::Residue *>::const_iterator it_set;
for (it_set=s.begin(); it_set!=s.end(); ++it_set) {
std::cout << " " << residue_spec_t(*it_set) << std::endl;
}
}
std::cout << "debug:: in init_from_residue_vec() here with these "
<< non_bonded_neighbour_residues.size() << " non-bonded neighbours"
<< std::endl;
for (std::size_t jj=0; jj<non_bonded_neighbour_residues.size(); jj++) {
std::cout << " " << residue_spec_t(non_bonded_neighbour_residues[jj]) << std::endl;
}
}
// std::cout << " DEBUG:: made " << bpc.size() << " bonded flanking pairs " << std::endl;
// passed and flanking
//
std::vector<mmdb::Residue *> all_residues;
std::vector<mmdb::Residue *>::const_iterator it;
for (unsigned int i=0; i<residues.size(); i++)
all_residues.push_back(residues[i].second);
// we don't need to calculate the NBC for all atoms - just these ones:
n_atoms_limit_for_nbc = 0;
for (unsigned int i=0; i<residues.size(); i++) {
n_atoms_limit_for_nbc += residues[i].second->GetNumberOfAtoms();
// std::cout << "Here in init_from_residue_vec() with n_atoms_limit_for_nbc " << n_atoms_limit_for_nbc << std::endl;
}
// Include only the fixed residues, because they are the flankers,
// the other residues are the ones in the passed residues vector.
// We don't have members of bonded_pair_container_t that are both
// fixed.
//
// 20151128 only include the residues once (test that they are not there first)
//
int n_bonded_flankers_in_total = 0; // debug/info counter
for (unsigned int i=0; i<bpc.size(); i++) {
if (bpc[i].is_fixed_first) {
it = std::find(all_residues.begin(),
all_residues.end(),
bpc[i].res_1);
if (it == all_residues.end()) {
all_residues.push_back(bpc[i].res_1);
n_bonded_flankers_in_total++;
}
}
if (bpc[i].is_fixed_second) {
it = std::find(all_residues.begin(),
all_residues.end(),
bpc[i].res_2);
if (it == all_residues.end()) {
all_residues.push_back(bpc[i].res_2);
n_bonded_flankers_in_total++;
}
}
}
// Finally add the neighbour residues that are not bonded:
for (unsigned int ires=0; ires<non_bonded_neighbour_residues.size(); ires++) {
it = std::find(all_residues.begin(),
all_residues.end(),
non_bonded_neighbour_residues[ires]);
if (it == all_residues.end())
all_residues.push_back(non_bonded_neighbour_residues[ires]);
}
if (debug) {
std::cout << " DEBUG:: There are " << residues.size() << " passed residues and "
<< all_residues.size() << " residues total (including flankers)"
<< " with " << non_bonded_neighbour_residues.size()
<< " non-bonded neighbours" << std::endl;
std::cout << "These are the non-bonded neighbour residues " << std::endl;
for (unsigned int i=0; i<non_bonded_neighbour_residues.size(); i++) {
std::cout << " " << residue_spec_t(non_bonded_neighbour_residues[i]) << std::endl;
}
}
n_atoms = 0;