1 //
2 // coor.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <set>
33
34 #include <util/misc/math.h>
35
36 #include <util/class/scexception.h>
37 #include <util/misc/formio.h>
38 #include <util/state/stateio.h>
39 #include <math/scmat/matrix.h>
40 #include <chemistry/molecule/molecule.h>
41 #include <chemistry/molecule/coor.h>
42 #include <chemistry/molecule/simple.h>
43 #include <chemistry/molecule/localdef.h>
44
45 #include <util/container/bitarray.h>
46
47 using namespace std;
48 using namespace sc;
49
50 ///////////////////////////////////////////////////////////////////////////
51 // members of IntCoor
52
53 double IntCoor::bohr_conv = 0.52917706;
54 double IntCoor::radian_conv = 180.0/M_PI;
55
56 static ClassDesc IntCoor_cd(
57 typeid(IntCoor),"IntCoor",1,"public SavableState",
58 0, 0, 0);
59
IntCoor(const char * re)60 IntCoor::IntCoor(const char *re):
61 label_(0), value_(0.0)
62 {
63 if (!re) re = "noname";
64 label_=new char[strlen(re)+1]; strcpy(label_,re);
65 }
66
IntCoor(const IntCoor & c)67 IntCoor::IntCoor(const IntCoor& c):
68 label_(0)
69 {
70 value_ = c.value_;
71 if (c.label_) label_ = strcpy(new char[strlen(c.label_)+1],c.label_);
72 }
73
IntCoor(const Ref<KeyVal> & keyval)74 IntCoor::IntCoor(const Ref<KeyVal>&keyval)
75 {
76 label_ = keyval->pcharvalue("label");
77 value_ = keyval->doublevalue("value");
78
79 if (keyval->exists("unit")) {
80 std::string unit(keyval->stringvalue("unit"));
81 if (unit == "bohr") {
82 }
83 else if (unit == "angstrom") {
84 value_ /= bohr_conv;
85 }
86 else if (unit == "radian") {
87 }
88 else if (unit == "degree") {
89 value_ *= M_PI/180.0;
90 }
91 else {
92 throw InputError("unrecognized unit value",
93 __FILE__, __LINE__,
94 "unit", unit.c_str(),
95 this->class_desc());
96 }
97 }
98 }
99
IntCoor(StateIn & si)100 IntCoor::IntCoor(StateIn& si):
101 SavableState(si)
102 {
103 si.get(value_);
104 si.getstring(label_);
105 }
106
~IntCoor()107 IntCoor::~IntCoor()
108 {
109 if (label_) delete[] label_;
110 }
111
112 void
save_data_state(StateOut & so)113 IntCoor::save_data_state(StateOut& so)
114 {
115 so.put(value_);
116 so.putstring(label_);
117 }
118
119 const char*
label() const120 IntCoor::label() const
121 {
122 return label_;
123 }
124
125 double
value() const126 IntCoor::value() const
127 {
128 return value_;
129 }
130
131 void
set_value(double v)132 IntCoor::set_value(double v)
133 {
134 value_ = v;
135 }
136
137 void
print(ostream & o) const138 IntCoor::print(ostream &o) const
139 {
140 print_details(0,o);
141 }
142
143 void
print_details(const Ref<Molecule> & mol,ostream & os) const144 IntCoor::print_details(const Ref<Molecule> &mol, ostream& os) const
145 {
146 os.setf(ios::fixed,ios::floatfield);
147 os.precision(10);
148 os.setf(ios::left,ios::adjustfield);
149 os.width(10);
150
151 os << indent
152 << scprintf("%-5s \"%10s\" %15.10f\n",ctype(),label(),preferred_value());
153 }
154
155 double
preferred_value() const156 IntCoor::preferred_value() const
157 {
158 return value_;
159 }
160
161 ///////////////////////////////////////////////////////////////////////////
162 // members of SetIntCoor
163
164 static ClassDesc SetIntCoor_cd(
165 typeid(SetIntCoor),"SetIntCoor",1,"public SavableState",
166 create<SetIntCoor>, create<SetIntCoor>, create<SetIntCoor>);
167
SetIntCoor()168 SetIntCoor::SetIntCoor()
169 {
170 }
171
SetIntCoor(const Ref<KeyVal> & keyval)172 SetIntCoor::SetIntCoor(const Ref<KeyVal>& keyval)
173 {
174 int n = keyval->count();
175
176 Ref<IntCoorGen> gen; gen << keyval->describedclassvalue("generator");
177
178 if (gen.null() && !n) {
179 throw InputError("not a vector and no generator given",
180 __FILE__, __LINE__,
181 0, 0,
182 class_desc());
183 }
184
185 if (gen.nonnull()) {
186 // Make sure that gen doesn't delete me before my reference
187 // count gets incremented.
188 this->reference();
189 gen->generate(this);
190 // Now it is safe to decrement my reference count back down to zero.
191 this->dereference();
192 }
193
194 for (int i=0; i<n; i++) {
195 Ref<IntCoor> coori; coori << keyval->describedclassvalue(i);
196 coor_.push_back(coori);
197 }
198 }
199
SetIntCoor(StateIn & s)200 SetIntCoor::SetIntCoor(StateIn& s):
201 SavableState(s)
202 {
203 int n;
204 s.get(n);
205
206 Ref<IntCoor> tmp;
207 for (int i=0; i<n; i++) {
208 tmp << SavableState::restore_state(s);
209 coor_.push_back(tmp);
210 }
211 }
212
~SetIntCoor()213 SetIntCoor::~SetIntCoor()
214 {
215 }
216
217 void
save_data_state(StateOut & s)218 SetIntCoor::save_data_state(StateOut& s)
219 {
220 int n = coor_.size();
221 s.put(n);
222
223 for (int i=0; i<n; i++) {
224 SavableState::save_state(coor_[i].pointer(),s);
225 }
226 }
227
228 void
add(const Ref<IntCoor> & coor)229 SetIntCoor::add(const Ref<IntCoor>& coor)
230 {
231 coor_.push_back(coor);
232 }
233
234 void
add(const Ref<SetIntCoor> & coor)235 SetIntCoor::add(const Ref<SetIntCoor>& coor)
236 {
237 for (int i=0; i<coor->n(); i++) {
238 coor_.push_back(coor->coor(i));
239 }
240 }
241
242 void
pop()243 SetIntCoor::pop()
244 {
245 coor_.pop_back();
246 }
247
248 int
n() const249 SetIntCoor::n() const
250 {
251 return coor_.size();
252 }
253
254 Ref<IntCoor>
coor(int i) const255 SetIntCoor::coor(int i) const
256 {
257 return coor_[i];
258 }
259
260 // compute the bmatrix by finite displacements
261 void
fd_bmat(const Ref<Molecule> & mol,RefSCMatrix & fd_bmatrix)262 SetIntCoor::fd_bmat(const Ref<Molecule>& mol,RefSCMatrix& fd_bmatrix)
263 {
264 Ref<SCMatrixKit> kit = fd_bmatrix.kit();
265
266 fd_bmatrix.assign(0.0);
267
268 int i;
269 Molecule& m = * mol.pointer();
270
271 const double cart_disp = 0.01;
272
273 RefSCDimension dn3(fd_bmatrix.coldim());
274 RefSCDimension dnc(fd_bmatrix.rowdim());
275 int n3 = dn3.n();
276 int nc = dnc.n();
277 RefSCVector internal(dnc,kit);
278 RefSCVector internal_p(dnc,kit);
279 RefSCVector internal_m(dnc,kit);
280
281 // the internal coordinates
282 update_values(mol);
283 for (i=0; i<nc; i++) {
284 internal(i) = coor_[i]->value();
285 }
286
287 // the finite displacement bmat
288 for (i=0; i<n3; i++) {
289 // the plus displacement
290 m.r(i/3,i%3) += cart_disp;
291 update_values(mol);
292 int j;
293 for (j=0; j<nc; j++) {
294 internal_p(j) = coor_[j]->value();
295 }
296 // the minus displacement
297 m.r(i/3,i%3) -= 2.0*cart_disp;
298 update_values(mol);
299 for (j=0; j<nc; j++) {
300 internal_m(j) = coor_[j]->value();
301 }
302 // reset the cartesian coordinate to its original value
303 m.r(i/3,i%3) += cart_disp;
304
305 // construct the entries in the finite displacement bmat
306 for (j=0; j<nc; j++) {
307 fd_bmatrix(j,i) = (internal_p(j)-internal_m(j))/(2.0*cart_disp);
308 }
309 }
310 }
311
312 void
bmat(const Ref<Molecule> & mol,RefSCMatrix & bmat)313 SetIntCoor::bmat(const Ref<Molecule>& mol, RefSCMatrix& bmat)
314 {
315 bmat.assign(0.0);
316
317 int i, ncoor = n();
318
319 RefSCVector bmatrow(bmat.coldim(),bmat.kit());
320 // send the rows of the b matrix to each of the coordinates
321 for (i=0; i<ncoor; i++) {
322 bmatrow.assign(0.0);
323 coor_[i]->bmat(mol,bmatrow);
324 bmat.assign_row(bmatrow,i);
325 }
326 }
327
328 void
guess_hessian(Ref<Molecule> & mol,RefSymmSCMatrix & hessian)329 SetIntCoor::guess_hessian(Ref<Molecule>& mol,RefSymmSCMatrix& hessian)
330 {
331 int ncoor = hessian.n();
332
333 hessian.assign(0.0);
334 for (int i=0; i<ncoor; i++) {
335 hessian(i,i) = coor_[i]->force_constant(mol);
336 }
337 }
338
339 void
print_details(const Ref<Molecule> & mol,ostream & os) const340 SetIntCoor::print_details(const Ref<Molecule> &mol, ostream& os) const
341 {
342 int i;
343
344 for(i=0; i<coor_.size(); i++) {
345 coor_[i]->print_details(mol,os);
346 }
347 }
348
349 void
update_values(const Ref<Molecule> & mol)350 SetIntCoor::update_values(const Ref<Molecule>&mol)
351 {
352 for (int i=0; i<coor_.size(); i++) {
353 coor_[i]->update_value(mol);
354 }
355 }
356
357 void
values_to_vector(const RefSCVector & v)358 SetIntCoor::values_to_vector(const RefSCVector&v)
359 {
360 for (int i=0; i<coor_.size(); i++) {
361 v(i) = coor_[i]->value();
362 }
363 }
364
365 void
clear()366 SetIntCoor::clear()
367 {
368 coor_.clear();
369 }
370
371 ///////////////////////////////////////////////////////////////////////////
372 // members of SumIntCoor
373
374 static ClassDesc SumIntCoor_cd(
375 typeid(SumIntCoor),"SumIntCoor",1,"public IntCoor",
376 0, create<SumIntCoor>, create<SumIntCoor>);
377
SumIntCoor(const char * label)378 SumIntCoor::SumIntCoor(const char* label):
379 IntCoor(label)
380 {
381 }
382
SumIntCoor(const Ref<KeyVal> & keyval)383 SumIntCoor::SumIntCoor(const Ref<KeyVal>&keyval):
384 IntCoor(keyval)
385 {
386 static const char* coor = "coor";
387 static const char* coef = "coef";
388 int n = keyval->count(coor);
389 int ncoef = keyval->count(coef);
390 if (n != ncoef) {
391 throw InputError("coor and coef do not have the same dimension",
392 __FILE__, __LINE__, 0, 0, class_desc());
393 }
394 if (!n) {
395 throw InputError("coor and coef are zero length",
396 __FILE__, __LINE__, 0, 0, class_desc());
397 }
398
399 for (int i=0; i<n; i++) {
400 double coe = keyval->doublevalue(coef,i);
401 Ref<IntCoor> coo; coo << keyval->describedclassvalue(coor,i);
402 add(coo,coe);
403 }
404 }
405
SumIntCoor(StateIn & s)406 SumIntCoor::SumIntCoor(StateIn&s):
407 IntCoor(s)
408 {
409 int n;
410 s.get(n);
411
412 coef_.resize(n);
413 coor_.resize(n);
414 for (int i=0; i<n; i++) {
415 s.get(coef_[i]);
416 coor_[i] << SavableState::restore_state(s);
417 }
418 }
419
~SumIntCoor()420 SumIntCoor::~SumIntCoor()
421 {
422 }
423
424 void
save_data_state(StateOut & s)425 SumIntCoor::save_data_state(StateOut&s)
426 {
427 int n = coef_.size();
428 IntCoor::save_data_state(s);
429 s.put(int(coef_.size()));
430
431 for (int i=0; i<n; i++) {
432 s.put(coef_[i]);
433 SavableState::save_state(coor_[i].pointer(),s);
434 }
435 }
436
437 int
n()438 SumIntCoor::n()
439 {
440 return coef_.size();
441 }
442
443 void
add(Ref<IntCoor> & coor,double coef)444 SumIntCoor::add(Ref<IntCoor>&coor,double coef)
445 {
446 // if a sum is added to a sum, unfold the nested sum
447 SumIntCoor* scoor = dynamic_cast<SumIntCoor*>(coor.pointer());
448 if (scoor) {
449 int l = scoor->coor_.size();
450 for (int i=0; i<l; i++) {
451 add(scoor->coor_[i],coef * scoor->coef_[i]);
452 }
453 }
454 else {
455 int l = coef_.size();
456 for (int i=0; i<l; i++) {
457 if (coor_[i]->equivalent(coor)) {
458 coef_[i] += coef;
459 return;
460 }
461 }
462 coef_.resize(l+1);
463 coor_.resize(l+1);
464 coef_[l] = coef;
465 coor_[l] = coor;
466 }
467 }
468
469 int
equivalent(Ref<IntCoor> & c)470 SumIntCoor::equivalent(Ref<IntCoor>&c)
471 {
472 return 0;
473 }
474
475 // this normalizes and makes the biggest coordinate positive
476 void
normalize()477 SumIntCoor::normalize()
478 {
479 int i;
480 int n = coef_.size();
481 double norm = 0.0;
482
483 double biggest = 0.0;
484 for (i=0; i<n; i++) {
485 norm += coef_[i] * coef_[i];
486 if (fabs(biggest) < fabs(coef_[i])) biggest = coef_[i];
487 }
488 norm = (biggest < 0.0? -1.0:1.0)/sqrt(norm);
489
490 for (i=0; i<n; i++) {
491 coef_[i] = coef_[i]*norm;
492 }
493 }
494
495 double
preferred_value() const496 SumIntCoor::preferred_value() const
497 {
498 return value_;
499 }
500
501 const char*
ctype() const502 SumIntCoor::ctype() const
503 {
504 return "SUM";
505 }
506
507 void
print_details(const Ref<Molecule> & mol,ostream & os) const508 SumIntCoor::print_details(const Ref<Molecule> &mol, ostream& os) const
509 {
510 int initial_indent = SCFormIO::getindent(os);
511 int i;
512
513 os << indent
514 << scprintf("%-5s %10s %14.10f\n",ctype(),
515 (label()?label():""), preferred_value());
516
517 for(i=0; i<coor_.size(); i++) {
518 os << incindent
519 << indent << scprintf("%14.10f ",coef_[i]);
520
521 SCFormIO::setindent(os, SCFormIO::getindent(os) + 15);
522 os << skipnextindent;
523 coor_[i]->print_details(mol,os);
524 SCFormIO::setindent(os, initial_indent);
525 }
526 }
527
528 // the SumIntCoor should be normalized before this is called.
529 double
force_constant(Ref<Molecule> & molecule)530 SumIntCoor::force_constant(Ref<Molecule>&molecule)
531 {
532 double fc = 0.0;
533
534 for (int i=0; i<n(); i++) {
535 fc += coef_[i] * coef_[i] * coor_[i]->force_constant(molecule);
536 }
537
538 return fc;
539 }
540
541 void
update_value(const Ref<Molecule> & molecule)542 SumIntCoor::update_value(const Ref<Molecule>&molecule)
543 {
544 int i, l = n();
545
546 value_ = 0.0;
547 for (i=0; i<l; i++) {
548 coor_[i]->update_value(molecule);
549 #if OLD_BMAT
550 if (dynamic_cast<StreSimpleCo*>(coor_[i]))
551 value_ += coef_[i] * dynamic_cast<StreSimpleCo*>(coor_[i])->angstrom();
552 else
553 #endif
554 value_ += coef_[i] * coor_[i]->value();
555 }
556 }
557
558 void
bmat(const Ref<Molecule> & molecule,RefSCVector & bmat,double coef)559 SumIntCoor::bmat(const Ref<Molecule>&molecule,RefSCVector&bmat,double coef)
560 {
561 int i, l = n();
562
563 for (i=0; i<l; i++) {
564 coor_[i]->bmat(molecule,bmat,coef*coef_[i]);
565 }
566 }
567
568 ///////////////////////////////////////////////////////////////////////////
569 // members of MolecularCoor
570
571 static ClassDesc MolecularCoor_cd(
572 typeid(MolecularCoor),"MolecularCoor",1,"public SavableState",
573 0, 0, 0);
574
MolecularCoor(Ref<Molecule> & mol)575 MolecularCoor::MolecularCoor(Ref<Molecule>&mol):
576 molecule_(mol)
577 {
578 debug_ = 0;
579 matrixkit_ = SCMatrixKit::default_matrixkit();
580 dnatom3_ = new SCDimension(3*molecule_->natom());
581 }
582
MolecularCoor(const Ref<KeyVal> & keyval)583 MolecularCoor::MolecularCoor(const Ref<KeyVal>&keyval)
584 {
585 molecule_ << keyval->describedclassvalue("molecule");
586
587 if (molecule_.null()) {
588 throw InputError("missing input", __FILE__, __LINE__,
589 "molecule", 0, class_desc());
590 }
591
592 debug_ = keyval->intvalue("debug");
593
594 matrixkit_ << keyval->describedclassvalue("matrixkit");
595 dnatom3_ << keyval->describedclassvalue("natom3");
596
597 if (matrixkit_.null()) matrixkit_ = SCMatrixKit::default_matrixkit();
598
599 if (dnatom3_.null()) dnatom3_ = new SCDimension(3*molecule_->natom());
600 else if (dnatom3_->n() != 3 * molecule_->natom()) {
601 throw InputError("natom3 given but not consistent with molecule",
602 __FILE__, __LINE__, "natom3", 0, class_desc());
603 }
604 }
605
MolecularCoor(StateIn & s)606 MolecularCoor::MolecularCoor(StateIn&s):
607 SavableState(s)
608 {
609 debug_ = 0;
610 matrixkit_ = SCMatrixKit::default_matrixkit();
611 molecule_ << SavableState::restore_state(s);
612 dnatom3_ << SavableState::restore_state(s);
613 }
614
~MolecularCoor()615 MolecularCoor::~MolecularCoor()
616 {
617 }
618
619 void
save_data_state(StateOut & s)620 MolecularCoor::save_data_state(StateOut&s)
621 {
622 SavableState::save_state(molecule_.pointer(),s);
623 SavableState::save_state(dnatom3_.pointer(),s);
624 }
625
626 int
nconstrained()627 MolecularCoor::nconstrained()
628 {
629 return 0;
630 }
631
632 // The default action is to never change the coordinates.
633 Ref<NonlinearTransform>
change_coordinates()634 MolecularCoor::change_coordinates()
635 {
636 return 0;
637 }
638
639 int
to_cartesian(const RefSCVector & internal)640 MolecularCoor::to_cartesian(const RefSCVector&internal)
641 {
642 return to_cartesian(molecule_, internal);
643 }
644
645 ///////////////////////////////////////////////////////////////////////////
646 // members of IntCoorGen
647
648 static ClassDesc IntCoorGen_cd(
649 typeid(IntCoorGen),"IntCoorGen",2,"public SavableState",
650 0, create<IntCoorGen>, create<IntCoorGen>);
651
IntCoorGen(const Ref<Molecule> & mol,int nextra_bonds,int * extra_bonds)652 IntCoorGen::IntCoorGen(const Ref<Molecule>& mol,
653 int nextra_bonds, int *extra_bonds)
654 {
655 init_constants();
656
657 molecule_ = mol;
658 nextra_bonds_ = nextra_bonds;
659 extra_bonds_ = extra_bonds;
660 }
661
IntCoorGen(const Ref<KeyVal> & keyval)662 IntCoorGen::IntCoorGen(const Ref<KeyVal>& keyval)
663 {
664 init_constants();
665
666 molecule_ << keyval->describedclassvalue("molecule");
667
668 radius_scale_factor_
669 = keyval->doublevalue("radius_scale_factor",
670 KeyValValuedouble(radius_scale_factor_));
671
672 // degrees
673 linear_bend_thres_
674 = keyval->doublevalue("linear_bend_threshold",
675 KeyValValuedouble(linear_bend_thres_));
676
677 // entered in degrees; stored as cos(theta)
678 linear_tors_thres_
679 = keyval->doublevalue("linear_tors_threshold",
680 KeyValValuedouble(linear_tors_thres_));
681
682 linear_bends_
683 = keyval->booleanvalue("linear_bend",
684 KeyValValueboolean(linear_bends_));
685
686 linear_lbends_
687 = keyval->booleanvalue("linear_lbend",
688 KeyValValueboolean(linear_lbends_));
689
690 linear_tors_
691 = keyval->booleanvalue("linear_tors",
692 KeyValValueboolean(linear_tors_));
693
694 linear_stors_
695 = keyval->booleanvalue("linear_stors",
696 KeyValValueboolean(linear_stors_));
697
698 // the extra_bonds list is given as a vector of atom numbers
699 // (atom numbering starts at 1)
700 nextra_bonds_ = keyval->count("extra_bonds");
701 nextra_bonds_ /= 2;
702 if (nextra_bonds_) {
703 extra_bonds_ = new int[nextra_bonds_*2];
704 for (int i=0; i<nextra_bonds_*2; i++) {
705 extra_bonds_[i] = keyval->intvalue("extra_bonds",i);
706 if (keyval->error() != KeyVal::OK) {
707 throw InputError("missing an expected integer value",
708 __FILE__, __LINE__, "extra_bonds", 0,
709 class_desc());
710 }
711 }
712 }
713 else {
714 extra_bonds_ = 0;
715 }
716 }
717
IntCoorGen(StateIn & s)718 IntCoorGen::IntCoorGen(StateIn& s):
719 SavableState(s)
720 {
721 molecule_ << SavableState::restore_state(s);
722 s.get(linear_bends_);
723 if (s.version(::class_desc<IntCoorGen>()) >= 2) {
724 s.get(linear_lbends_);
725 }
726 s.get(linear_tors_);
727 s.get(linear_stors_);
728 s.get(linear_bend_thres_);
729 s.get(linear_tors_thres_);
730 s.get(nextra_bonds_);
731 s.get(extra_bonds_);
732 s.get(radius_scale_factor_);
733 }
734
735 void
init_constants()736 IntCoorGen::init_constants()
737 {
738 nextra_bonds_ = 0;
739 extra_bonds_ = 0;
740 radius_scale_factor_ = 1.1;
741 linear_bend_thres_ = 1.0;
742 linear_tors_thres_ = 1.0;
743 linear_bends_ = 0;
744 linear_lbends_ = 1;
745 linear_tors_ = 0;
746 linear_stors_ = 1;
747 }
748
~IntCoorGen()749 IntCoorGen::~IntCoorGen()
750 {
751 if (extra_bonds_) delete[] extra_bonds_;
752 }
753
754 void
save_data_state(StateOut & s)755 IntCoorGen::save_data_state(StateOut& s)
756 {
757 SavableState::save_state(molecule_.pointer(),s);
758 s.put(linear_bends_);
759 s.put(linear_lbends_);
760 s.put(linear_tors_);
761 s.put(linear_stors_);
762 s.put(linear_bend_thres_);
763 s.put(linear_tors_thres_);
764 s.put(nextra_bonds_);
765 s.put(extra_bonds_,2*nextra_bonds_);
766 s.put(radius_scale_factor_);
767 }
768
769 void
print(ostream & out) const770 IntCoorGen::print(ostream& out) const
771 {
772 out << indent << "IntCoorGen:" << endl << incindent
773 << indent << "linear_bends = " << linear_bends_ << endl
774 << indent << "linear_lbends = " << linear_lbends_ << endl
775 << indent << "linear_tors = " << linear_tors_ << endl
776 << indent << "linear_stors = " << linear_stors_ << endl
777 << indent << scprintf("linear_bend_threshold = %f\n",linear_bend_thres_)
778 << indent << scprintf("linear_tors_threshold = %f\n",linear_tors_thres_)
779 << indent << scprintf("radius_scale_factor = %f\n",radius_scale_factor_)
780 << indent << "nextra_bonds = " << nextra_bonds_ << endl
781 << decindent;
782 }
783
784 static void
find_bonds(Molecule & m,BitArrayLTri & bonds,double radius_scale_factor_)785 find_bonds(Molecule &m, BitArrayLTri &bonds,
786 double radius_scale_factor_)
787 {
788 int i, j;
789 for(i=0; i < m.natom(); i++) {
790 double at_rad_i = m.atominfo()->atomic_radius(m.Z(i));
791 SCVector3 ri(m.r(i));
792
793 for(j=0; j < i; j++) {
794 double at_rad_j = m.atominfo()->atomic_radius(m.Z(j));
795 SCVector3 rj(m.r(j));
796
797 if (ri.dist(rj)
798 < radius_scale_factor_*(at_rad_i+at_rad_j))
799 bonds.set(i,j);
800 }
801 }
802
803 // check for groups of atoms bound to nothing
804 std::set<int> boundatoms;
805 std::set<int> newatoms, nextnewatoms;
806 // start out with atom 0
807 newatoms.insert(0);
808 std::set<int>::iterator iatom;
809 int warning_printed = 0;
810 while (newatoms.size() > 0) {
811 while (newatoms.size() > 0) {
812 // newatoms gets merged into boundatoms
813 for (iatom=newatoms.begin(); iatom!=newatoms.end(); iatom++) {
814 boundatoms.insert(*iatom);
815 }
816 // set nextnewatoms to atoms bound to boundatoms that are not already
817 // in boundatoms
818 nextnewatoms.clear();
819 for (iatom=newatoms.begin(); iatom!=newatoms.end(); iatom++) {
820 int atom = *iatom;
821 for (i=0; i<m.natom(); i++) {
822 if (bonds(i,atom) && boundatoms.find(i) == boundatoms.end()) {
823 nextnewatoms.insert(i);
824 }
825 }
826 }
827 // set newatoms to nextnewatoms to start off the next iteration
828 newatoms.clear();
829 for (iatom=nextnewatoms.begin(); iatom!=nextnewatoms.end(); iatom++) {
830 newatoms.insert(*iatom);
831 }
832 }
833
834 if (boundatoms.size() != m.natom()) {
835 if (!warning_printed) {
836 warning_printed = 1;
837 ExEnv::out0()
838 << indent << "WARNING: two unbound groups of atoms" << endl
839 << indent << " consider using extra_bonds input" << endl
840 << endl;
841 }
842 // find an unbound group
843 double nearest_dist;
844 int nearest_bound = -1, nearest_unbound = -1;
845 for(i=0; i < m.natom(); i++) {
846 if (boundatoms.find(i) == boundatoms.end()) {
847 SCVector3 ri(m.r(i));
848 for (iatom=boundatoms.begin(); iatom!=boundatoms.end(); iatom++) {
849 SCVector3 rj(m.r(*iatom));
850 double d = ri.dist(rj);
851 if (nearest_unbound == -1 || d < nearest_dist) {
852 nearest_dist = d;
853 nearest_bound = *iatom;
854 nearest_unbound = i;
855 }
856 }
857 }
858 }
859 if (nearest_bound == -1) {
860 throw ProgrammingError("impossible error generating coordinates",
861 __FILE__, __LINE__);
862 }
863 // add all bound atoms within a certain distance of nearest_unbound
864 // --- should really do this for all atoms that nearest_unbound
865 // --- is already bound to, too
866 SCVector3 rnearest_unbound(m.r(nearest_unbound));
867 for (iatom=boundatoms.begin(); iatom!=boundatoms.end(); iatom++) {
868 SCVector3 r(m.r(*iatom));
869 if (*iatom == nearest_bound
870 || rnearest_unbound.dist(r) < 1.1 * nearest_dist) {
871 ExEnv::out0() << indent
872 << " adding bond between "
873 << *iatom+1 << " and " << nearest_unbound+1 << endl;
874 bonds.set(*iatom,nearest_unbound);
875 }
876 }
877 newatoms.insert(nearest_unbound);
878 }
879 }
880 }
881
882 void
generate(const Ref<SetIntCoor> & sic)883 IntCoorGen::generate(const Ref<SetIntCoor>& sic)
884 {
885 int i;
886 Molecule& m = *molecule_.pointer();
887
888 // let's go through the geometry and find all the close contacts
889 // bonds is a lower triangle matrix of 1's and 0's indicating whether
890 // there is a bond between atoms i and j
891
892 BitArrayLTri bonds(m.natom(),m.natom());
893
894 for (i=0; i<nextra_bonds_; i++) {
895 bonds.set(extra_bonds_[i*2]-1,extra_bonds_[i*2+1]-1);
896 }
897
898 find_bonds(m, bonds, radius_scale_factor_);
899
900 // compute the simple internal coordinates by type
901 add_bonds(sic,bonds,m);
902 add_bends(sic,bonds,m);
903 add_tors(sic,bonds,m);
904 add_out(sic,bonds,m);
905
906 ExEnv::out0() << endl << indent
907 << "IntCoorGen: generated " << sic->n() << " coordinates." << endl;
908 }
909
910 ///////////////////////////////////////////////////////////////////////////
911 // auxillary functions of IntCoorGen
912
913 /*
914 * the following are translations of functions written by Gregory Humphreys
915 * at the NIH
916 */
917
918 /*
919 * for each bonded pair, add an entry to the simple coord list
920 */
921
922 void
add_bonds(const Ref<SetIntCoor> & list,BitArrayLTri & bonds,Molecule & m)923 IntCoorGen::add_bonds(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m)
924 {
925 int i,j,ij;
926 int labelc=0;
927 char label[80];
928
929 for(i=ij=0; i < m.natom(); i++) {
930 for(j=0; j <= i; j++,ij++) {
931 if(bonds[ij]) {
932 labelc++;
933 sprintf(label,"s%d",labelc);
934 list->add(new Stre(label,j+1,i+1));
935 }
936 }
937 }
938 }
939
940 /*
941 * return 1 if all three atoms are nearly on the same line.
942 */
943
944 // returns fabs(cos(theta_ijk))
945 double
cos_ijk(Molecule & m,int i,int j,int k)946 IntCoorGen::cos_ijk(Molecule& m, int i, int j, int k)
947 {
948 SCVector3 a, b, c;
949 int xyz;
950 for (xyz=0; xyz<3; xyz++) {
951 a[xyz] = m.r(i,xyz);
952 b[xyz] = m.r(j,xyz);
953 c[xyz] = m.r(k,xyz);
954 }
955 SCVector3 ab = a - b;
956 SCVector3 cb = c - b;
957 return fabs(ab.dot(cb)/(ab.norm()*cb.norm()));
958 }
959
960 void
add_bends(const Ref<SetIntCoor> & list,BitArrayLTri & bonds,Molecule & m)961 IntCoorGen::add_bends(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m)
962 {
963 int i,j,k;
964 int labelc=0;
965 char label[80];
966
967 int n = m.natom();
968
969 double thres = cos(linear_bend_thres_*M_PI/180.0);
970
971 for(i=0; i < n; i++) {
972 SCVector3 ri(m.r(i));
973 for(j=0; j < n; j++) {
974 if(bonds(i,j)) {
975 SCVector3 rj(m.r(j));
976 for(k=0; k < i; k++) {
977 if(bonds(j,k)) {
978 SCVector3 rk(m.r(k));
979 int is_linear = (cos_ijk(m,i,j,k) >= thres);
980 if (linear_bends_ || !is_linear) {
981 labelc++;
982 sprintf(label,"b%d",labelc);
983 list->add(new Bend(label,k+1,j+1,i+1));
984 }
985 if (linear_lbends_ && is_linear) {
986 // find a unit vector roughly perp to the bonds
987 SCVector3 u;
988 // first try to find another atom, that'll help keep one of
989 // the coordinates totally symmetric in some cases
990 int most_perp_atom = -1;
991 double cos_most_perp = thres;
992 for (int l=0; l < n; l++) {
993 if (l == i || l == j || l == k) continue;
994 double tmp = cos_ijk(m,i,j,l);
995 if (tmp < cos_most_perp) {
996 cos_most_perp = tmp;
997 most_perp_atom = l;
998 }
999 }
1000 if (most_perp_atom != -1) {
1001 SCVector3 rmpa(m.r(most_perp_atom));
1002 u = rj-rmpa;
1003 u.normalize();
1004 }
1005 else {
1006 SCVector3 b1, b2;
1007 b1 = ri-rj;
1008 b2 = rk-rj;
1009 u = b1.perp_unit(b2);
1010 }
1011 labelc++;
1012 sprintf(label,"b%d",labelc);
1013 list->add(new LinIP(label,k+1,j+1,i+1,u));
1014 labelc++;
1015 sprintf(label,"b%d",labelc);
1016 list->add(new LinOP(label,k+1,j+1,i+1,u));
1017 }
1018 }
1019 }
1020 }
1021 }
1022 }
1023 }
1024
1025 /*
1026 * for each pair of bends which share a common bond, add a torsion
1027 */
1028
1029 /*
1030 * just look at the heavy-atom skeleton. return true if i is a terminal
1031 * atom.
1032 */
1033
1034 int
hterminal(Molecule & m,BitArrayLTri & bonds,int i)1035 IntCoorGen::hterminal(Molecule& m, BitArrayLTri& bonds, int i)
1036 {
1037 int nh=0;
1038 for (int j=0; j < m.natom(); j++)
1039 if (bonds(i,j) && m.Z(j) > 1) nh++;
1040 return (nh==1);
1041 }
1042
1043 void
add_tors(const Ref<SetIntCoor> & list,BitArrayLTri & bonds,Molecule & m)1044 IntCoorGen::add_tors(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m)
1045 {
1046 int i,j,k,l;
1047 int labelc=0;
1048 char label[80];
1049
1050 int n = m.natom();
1051
1052 double thres = cos(linear_tors_thres_*M_PI/180.0);
1053
1054 for(j=0; j < n; j++) {
1055 for(k=0; k < j; k++) {
1056 if(bonds(j,k)) {
1057 for(i=0; i < n; i++) {
1058 if(k==i) continue;
1059
1060 // no hydrogen torsions, ok?
1061 if (m.Z(i) == 1 && !hterminal(m,bonds,j)) continue;
1062
1063 if (bonds(j,i)) {
1064 int is_linear = 0;
1065 if (cos_ijk(m,i,j,k)>=thres) is_linear = 1;
1066
1067 for (l=0; l < n; l++) {
1068 if (l==j || l==i) continue;
1069
1070 // no hydrogen torsions, ok?
1071 if (m.Z(l) == 1 && !hterminal(m,bonds,k))
1072 continue;
1073
1074 if (bonds(k,l)) {
1075 if(cos_ijk(m,j,k,l)>=thres) is_linear = 1;
1076
1077 if (is_linear && linear_stors_) {
1078 labelc++;
1079 sprintf(label,"st%d",labelc);
1080 list->add(new ScaledTors(label,l+1,k+1,j+1,i+1));
1081 }
1082 if (!is_linear || linear_tors_) {
1083 labelc++;
1084 sprintf(label,"t%d",labelc);
1085 list->add(new Tors(label,l+1,k+1,j+1,i+1));
1086 }
1087 }
1088 }
1089 }
1090 }
1091 }
1092 }
1093 }
1094 }
1095
1096 void
add_out(const Ref<SetIntCoor> & list,BitArrayLTri & bonds,Molecule & m)1097 IntCoorGen::add_out(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m)
1098 {
1099 int i,j,k,l;
1100 int labelc=0;
1101 char label[80];
1102
1103 int n = m.natom();
1104
1105 // first find all tri-coordinate atoms
1106 for(i=0; i < n; i++) {
1107 if(bonds.degree(i)!=3) continue;
1108
1109 // then look for terminal atoms connected to i
1110 for(j=0; j < n; j++) {
1111 if(bonds(i,j) && bonds.degree(j)==1) {
1112
1113 for(k=0; k < n; k++) {
1114 if(k!=j && bonds(i,k)) {
1115 for(l=0; l < k; l++) {
1116 if(l!=j && bonds(i,l)) {
1117 labelc++;
1118 sprintf(label,"o%d",labelc);
1119 list->add(new Out(label,j+1,i+1,l+1,k+1));
1120 }
1121 }
1122 }
1123 }
1124 }
1125 }
1126 }
1127 }
1128
1129 int
nearest_contact(int i,Molecule & m)1130 IntCoorGen::nearest_contact(int i, Molecule& m)
1131 {
1132 double d=-1.0;
1133 int n=0;
1134
1135 SCVector3 ri(m.r(i));
1136
1137 for (int j=0; j < m.natom(); j++) {
1138 SCVector3 rj(m.r(j));
1139 double td = ri.dist(rj);
1140 if (j==i)
1141 continue;
1142 else if (d < 0 || td < d) {
1143 d = td;
1144 n = j;
1145 }
1146 }
1147
1148 return n;
1149 }
1150
1151 /////////////////////////////////////////////////////////////////////////////
1152
1153 // Local Variables:
1154 // mode: c++
1155 // c-file-style: "CLJ-CONDENSED"
1156 // End:
1157