1 //
2 // molecule.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 <math.h>
33 #include <string.h>
34 #include <stdio.h>
35
36 #include <util/class/scexception.h>
37 #include <util/misc/formio.h>
38 #include <util/state/stateio.h>
39 #include <chemistry/molecule/molecule.h>
40 #include <chemistry/molecule/formula.h>
41 #include <chemistry/molecule/localdef.h>
42 #include <math/scmat/cmatrix.h>
43
44 using namespace std;
45 using namespace sc;
46
47 //////////////////////////////////////////////////////////////////////////
48 // Molecule
49
50 static ClassDesc Molecule_cd(
51 typeid(Molecule),"Molecule",6,"public SavableState",
52 create<Molecule>, create<Molecule>, create<Molecule>);
53
Molecule()54 Molecule::Molecule():
55 natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
56 {
57 pg_ = new PointGroup;
58 atominfo_ = new AtomInfo();
59 geometry_units_ = new Units("bohr");
60 nuniq_ = 0;
61 equiv_ = 0;
62 nequiv_ = 0;
63 atom_to_uniq_ = 0;
64 q_Z_ = atominfo_->string_to_Z("Q");
65 include_q_ = false;
66 include_qq_ = false;
67 init_symmetry_info();
68 }
69
Molecule(const Molecule & mol)70 Molecule::Molecule(const Molecule& mol):
71 natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
72 {
73 nuniq_ = 0;
74 equiv_ = 0;
75 nequiv_ = 0;
76 atom_to_uniq_ = 0;
77 *this=mol;
78 }
79
~Molecule()80 Molecule::~Molecule()
81 {
82 clear();
83 }
84
85 void
clear()86 Molecule::clear()
87 {
88 if (r_) {
89 delete[] r_[0];
90 delete[] r_;
91 r_ = 0;
92 }
93 if (labels_) {
94 for (int i=0; i<natoms_; i++) {
95 delete[] labels_[i];
96 }
97 delete[] labels_;
98 labels_ = 0;
99 }
100 delete[] charges_;
101 charges_ = 0;
102 delete[] mass_;
103 mass_ = 0;
104 delete[] Z_;
105 Z_ = 0;
106
107 clear_symmetry_info();
108 }
109
110 void
throw_if_atom_duplicated(int begin,double tol)111 Molecule::throw_if_atom_duplicated(int begin, double tol)
112 {
113 for (int i=begin; i<natoms_; i++) {
114 SCVector3 ri(r_[i]);
115 for (int j=0; j<i; j++) {
116 SCVector3 rj(r_[j]);
117 if (ri.dist(rj) < tol) {
118 throw InputError("duplicated atom coordinate",
119 __FILE__, __LINE__, 0, 0, class_desc());
120 }
121 }
122 }
123 }
124
Molecule(const Ref<KeyVal> & input)125 Molecule::Molecule(const Ref<KeyVal>&input):
126 natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
127 {
128 nuniq_ = 0;
129 equiv_ = 0;
130 nequiv_ = 0;
131 atom_to_uniq_ = 0;
132
133 KeyValValueboolean kvfalse(0);
134 include_q_ = input->booleanvalue("include_q",kvfalse);
135 include_qq_ = input->booleanvalue("include_qq",kvfalse);
136
137 atominfo_ << input->describedclassvalue("atominfo");
138 if (atominfo_.null()) atominfo_ = new AtomInfo;
139 q_Z_ = atominfo_->string_to_Z("Q");
140 if (input->exists("pdb_file")) {
141 geometry_units_ = new Units("angstrom");
142 char* filename = input->pcharvalue("pdb_file");
143 read_pdb(filename);
144 delete[] filename;
145 }
146 else {
147 // check for old style units input first
148 if (input->booleanvalue("angstrom")
149 ||input->booleanvalue("angstroms")
150 ||input->booleanvalue("aangstrom")
151 ||input->booleanvalue("aangstroms")) {
152 geometry_units_ = new Units("angstrom");
153 }
154 // check for new style units input
155 else {
156 geometry_units_ = new Units(input->pcharvalue("unit"),
157 Units::Steal);
158 }
159
160 double conv = geometry_units_->to_atomic_units();
161
162 // get the number of atoms and make sure that the geometry and the
163 // atoms array have the same number of atoms.
164 // right now we read in the unique atoms...then we will symmetrize.
165 // the length of atoms must still equal the length of geometry, but
166 // we'll try to set up atom_labels such that different lengths are
167 // possible
168 int natom = input->count("geometry");
169 if (natom != input->count("atoms")) {
170 throw InputError("size of \"geometry\" != size of \"atoms\"",
171 __FILE__, __LINE__, 0, 0, class_desc());
172 }
173
174 int i;
175 for (i=0; i<natom; i++) {
176 char *name, *label;
177 int ghost = input->booleanvalue("ghost",i);
178 double charge = input->doublevalue("charge",i);
179 int have_charge = input->error() == KeyVal::OK;
180 if (ghost) {
181 have_charge = 1;
182 charge = 0.0;
183 }
184 add_atom(atominfo_->string_to_Z(name = input->pcharvalue("atoms",i)),
185 input->doublevalue("geometry",i,0)*conv,
186 input->doublevalue("geometry",i,1)*conv,
187 input->doublevalue("geometry",i,2)*conv,
188 label = input->pcharvalue("atom_labels",i),
189 input->doublevalue("mass",i),
190 have_charge, charge
191 );
192 delete[] name;
193 delete[] label;
194 }
195 }
196
197 char *symmetry = input->pcharvalue("symmetry");
198 double symtol = input->doublevalue("symmetry_tolerance",
199 KeyValValuedouble(1.0e-4));
200 nuniq_ = 0;
201 equiv_ = 0;
202 nequiv_ = 0;
203 atom_to_uniq_ = 0;
204 if (symmetry && !strcmp(symmetry,"auto")) {
205 set_point_group(highest_point_group(symtol), symtol*10.0);
206 }
207 else {
208 pg_ = new PointGroup(input);
209
210 // translate to the origin of the symmetry frame
211 double r[3];
212 for (int i=0; i<3; i++) {
213 r[i] = -pg_->origin()[i];
214 pg_->origin()[i] = 0;
215 }
216 translate(r);
217
218 if (input->booleanvalue("redundant_atoms")) {
219 init_symmetry_info();
220 cleanup_molecule(symtol);
221 }
222 else {
223 symmetrize();
224 // In case we were given redundant atoms, clean up
225 // the geometry so the symmetry is exact.
226 cleanup_molecule(symtol);
227 }
228 }
229 delete[] symmetry;
230 }
231
232 Molecule&
operator =(const Molecule & mol)233 Molecule::operator=(const Molecule& mol)
234 {
235 clear();
236
237 pg_ = new PointGroup(*(mol.pg_.pointer()));
238 atominfo_ = mol.atominfo_;
239 geometry_units_ = new Units(mol.geometry_units_->string_rep());
240
241 q_Z_ = mol.q_Z_;
242 include_q_ = mol.include_q_;
243 include_qq_ = mol.include_qq_;
244 q_atoms_ = mol.q_atoms_;
245 non_q_atoms_ = mol.non_q_atoms_;
246
247 natoms_ = mol.natoms_;
248
249 if (natoms_) {
250 if (mol.mass_) {
251 mass_ = new double[natoms_];
252 memcpy(mass_,mol.mass_,natoms_*sizeof(double));
253 }
254 if (mol.charges_) {
255 charges_ = new double[natoms_];
256 memcpy(charges_,mol.charges_,natoms_*sizeof(double));
257 }
258 if (mol.labels_) {
259 labels_ = new char *[natoms_];
260 for (int i=0; i<natoms_; i++) {
261 if (mol.labels_[i]) {
262 labels_[i] = strcpy(new char[strlen(mol.labels_[i])+1],
263 mol.labels_[i]);
264 }
265 else labels_[i] = 0;
266 }
267 }
268 r_ = new double*[natoms_];
269 r_[0] = new double[natoms_*3];
270 for (int i=0; i<natoms_; i++) {
271 r_[i] = &(r_[0][i*3]);
272 }
273 memcpy(r_[0], mol.r_[0], natoms_*3*sizeof(double));
274 Z_ = new int[natoms_];
275 memcpy(Z_, mol.Z_, natoms_*sizeof(int));
276 }
277
278 init_symmetry_info();
279
280 return *this;
281 }
282
283 void
add_atom(int Z,double x,double y,double z,const char * label,double mass,int have_charge,double charge)284 Molecule::add_atom(int Z,double x,double y,double z,
285 const char *label,double mass,
286 int have_charge, double charge)
287 {
288 int i;
289
290 // allocate new arrays
291 int *newZ = new int[natoms_+1];
292 double **newr = new double*[natoms_+1];
293 double *newr0 = new double[(natoms_+1)*3];
294 char **newlabels = 0;
295 if (label || labels_) {
296 newlabels = new char*[natoms_+1];
297 }
298 double *newcharges = 0;
299 if (have_charge || charges_) {
300 newcharges = new double[natoms_+1];
301 }
302 double *newmass = 0;
303 if (mass_ || mass != 0.0) {
304 newmass = new double[natoms_+1];
305 }
306
307 // setup the r_ pointers
308 for (i=0; i<=natoms_; i++) {
309 newr[i] = &(newr0[i*3]);
310 }
311
312 // copy old data to new arrays
313 if (natoms_) {
314 memcpy(newZ,Z_,sizeof(int)*natoms_);
315 memcpy(newr0,r_[0],sizeof(double)*natoms_*3);
316 if (labels_) {
317 memcpy(newlabels,labels_,sizeof(char*)*natoms_);
318 }
319 else if (newlabels) {
320 memset(newlabels,0,sizeof(char*)*natoms_);
321 }
322 if (charges_) {
323 memcpy(newcharges,charges_,sizeof(double)*natoms_);
324 }
325 else if (newcharges) {
326 for (i=0; i<natoms_; i++) newcharges[i] = Z_[i];
327 }
328 if (mass_) {
329 memcpy(newmass,mass_,sizeof(double)*natoms_);
330 }
331 else if (newmass) {
332 memset(newmass,0,sizeof(double)*natoms_);
333 }
334 }
335
336 // delete old data
337 delete[] Z_;
338 if (r_) {
339 delete[] r_[0];
340 delete[] r_;
341 }
342 delete[] labels_;
343 delete[] charges_;
344 delete[] mass_;
345
346 // setup new pointers
347 Z_ = newZ;
348 r_ = newr;
349 labels_ = newlabels;
350 charges_ = newcharges;
351 mass_ = newmass;
352
353 // copy info for this atom into arrays
354 Z_[natoms_] = Z;
355 r_[natoms_][0] = x;
356 r_[natoms_][1] = y;
357 r_[natoms_][2] = z;
358 if (mass_) mass_[natoms_] = mass;
359 if (label) {
360 labels_[natoms_] = strcpy(new char[strlen(label)+1],label);
361 }
362 else if (labels_) {
363 labels_[natoms_] = 0;
364 }
365 if (have_charge) {
366 charges_[natoms_] = charge;
367 }
368 else if (charges_) {
369 charges_[natoms_] = Z;
370 }
371
372 if (Z == q_Z_) {
373 q_atoms_.push_back(natoms_);
374 }
375 else {
376 non_q_atoms_.push_back(natoms_);
377 }
378
379 natoms_++;
380
381 throw_if_atom_duplicated(natoms_-1);
382 }
383
384 void
print_parsedkeyval(ostream & os,int print_pg,int print_unit,int number_atoms) const385 Molecule::print_parsedkeyval(ostream& os,
386 int print_pg,
387 int print_unit,
388 int number_atoms) const
389 {
390 int i;
391
392 double conv = geometry_units_->from_atomic_units();
393
394 if (print_pg) pg_->print(os);
395 if (print_unit && geometry_units_->string_rep()) {
396 os << indent
397 << "unit = \"" << geometry_units_->string_rep() << "\""
398 << endl;
399 }
400 os << indent << "{";
401 if (number_atoms) os << scprintf("%3s", "n");
402 os << scprintf(" %5s", "atoms");
403 if (labels_) os << scprintf(" %11s", "atom_labels");
404 int int_charges = 1;
405 if (charges_) {
406 for (i=0;i<natom();i++) if (charges_[i]!=(int)charges_[i]) int_charges=0;
407 if (int_charges) {
408 os << scprintf(" %7s", "charge");
409 }
410 else {
411 os << scprintf(" %17s", "charge");
412 }
413 }
414 os << scprintf(" %16s", "")
415 << scprintf(" %16s", "geometry ")
416 << scprintf(" %16s ", "");
417 os << "}={" << endl;
418 for (i=0; i<natom(); i++) {
419 os << indent;
420 if (number_atoms) os << scprintf(" %3d", i+1);
421 std::string symbol(atom_symbol(i));
422 os << scprintf(" %5s", symbol.c_str());
423 if (labels_) {
424 const char *lab = labels_[i];
425 if (lab == 0) lab = "";
426 char *qlab = new char[strlen(lab)+3];
427 strcpy(qlab,"\"");
428 strcat(qlab,lab);
429 strcat(qlab,"\"");
430 os << scprintf(" %11s",qlab);
431 delete[] qlab;
432 }
433 if (charges_) {
434 if (int_charges) os << scprintf(" %7.4f", charges_[i]);
435 else os << scprintf(" %17.15f", charges_[i]);
436 }
437 os << scprintf(" [% 16.10f", conv * r(i,0))
438 << scprintf(" % 16.10f", conv * r(i,1))
439 << scprintf(" % 16.10f]", conv * r(i,2))
440 << endl;
441 }
442 os << indent << "}" << endl;
443 }
444
445 void
print(ostream & os) const446 Molecule::print(ostream& os) const
447 {
448 int i;
449
450 MolecularFormula *mf = new MolecularFormula(this);
451 os << indent
452 << "Molecular formula: " << mf->formula() << endl;
453 delete mf;
454
455 os << indent << "molecule<Molecule>: (" << endl;
456 os << incindent;
457 print_parsedkeyval(os);
458 os << decindent;
459 os << indent << ")" << endl;
460
461 os << indent << "Atomic Masses:" << endl;
462 for (i=0; i<natom(); i+=5) {
463 os << indent;
464 for (int j=i; j<i+5 && j<natom(); j++) {
465 os << scprintf(" %10.5f", mass(j));
466 }
467 os << endl;
468 }
469 }
470
471 int
atom_label_to_index(const char * l) const472 Molecule::atom_label_to_index(const char *l) const
473 {
474 int i;
475 for (i=0; i<natom(); i++) {
476 if (label(i) && !strcmp(l,label(i))) return i;
477 }
478 return -1;
479 }
480
481 double*
charges() const482 Molecule::charges() const
483 {
484 double*result = new double[natoms_];
485 if (charges_) {
486 memcpy(result, charges_, sizeof(double)*natom());
487 }
488 else {
489 for (int i=0; i<natom(); i++) result[i] = Z_[i];
490 }
491 return result;
492 }
493
494 double
charge(int iatom) const495 Molecule::charge(int iatom) const
496 {
497 if (charges_) return charges_[iatom];
498 return Z_[iatom];
499 }
500
501 double
nuclear_charge() const502 Molecule::nuclear_charge() const
503 {
504 double c = 0.0;
505 if (include_q_) {
506 for (int i=0; i<natom(); i++) {
507 c += charge(i);
508 }
509 }
510 else {
511 for (int ii=0; ii<non_q_atoms_.size(); ii++) {
512 c += charge(non_q_atoms_[ii]);
513 }
514 }
515 return c;
516 }
517
save_data_state(StateOut & so)518 void Molecule::save_data_state(StateOut& so)
519 {
520 so.put(include_q_);
521 so.put(include_qq_);
522 so.put(natoms_);
523 SavableState::save_state(pg_.pointer(),so);
524 SavableState::save_state(geometry_units_.pointer(),so);
525 SavableState::save_state(atominfo_.pointer(),so);
526 if (natoms_) {
527 so.put(Z_, natoms_);
528 so.put_array_double(r_[0], natoms_*3);
529 so.put(charges_,natoms_);
530 }
531 if (mass_) {
532 so.put(1);
533 so.put_array_double(mass_, natoms_);
534 }
535 else {
536 so.put(0);
537 }
538 if (labels_){
539 so.put(1);
540 for (int i=0; i<natoms_; i++) {
541 so.putstring(labels_[i]);
542 }
543 }
544 else {
545 so.put(0);
546 }
547 }
548
Molecule(StateIn & si)549 Molecule::Molecule(StateIn& si):
550 SavableState(si),
551 natoms_(0), r_(0), Z_(0), mass_(0), labels_(0)
552 {
553 if (si.version(::class_desc<Molecule>()) < 4) {
554 throw FileOperationFailed("cannot restore from old molecules",
555 __FILE__, __LINE__, 0,
556 FileOperationFailed::Corrupt,
557 class_desc());
558 }
559 if (si.version(::class_desc<Molecule>()) < 6) {
560 include_q_ = false;
561 include_qq_ = false;
562 }
563 else {
564 si.get(include_q_);
565 si.get(include_qq_);
566 }
567 si.get(natoms_);
568 pg_ << SavableState::restore_state(si);
569 geometry_units_ << SavableState::restore_state(si);
570 atominfo_ << SavableState::restore_state(si);
571 if (natoms_) {
572 si.get(Z_);
573 r_ = new double*[natoms_];
574 r_[0] = new double[natoms_*3];
575 si.get_array_double(r_[0],natoms_*3);
576 for (int i=1; i<natoms_; i++) {
577 r_[i] = &(r_[0][i*3]);
578 }
579 if (si.version(::class_desc<Molecule>()) > 4) {
580 si.get(charges_);
581 }
582 else {
583 charges_ = 0;
584 }
585 }
586 int test;
587 si.get(test);
588 if (test) {
589 mass_ = new double[natoms_];
590 si.get_array_double(mass_, natoms_);
591 }
592 si.get(test);
593 if (test){
594 labels_ = new char*[natoms_];
595 for (int i=0; i<natoms_; i++) {
596 si.getstring(labels_[i]);
597 }
598 }
599
600 for (int i=0; i<natoms_; i++) {
601 if (Z_[i] == q_Z_) {
602 q_atoms_.push_back(i);
603 }
604 else {
605 non_q_atoms_.push_back(i);
606 }
607 }
608
609 nuniq_ = 0;
610 equiv_ = 0;
611 nequiv_ = 0;
612 atom_to_uniq_ = 0;
613 init_symmetry_info();
614 }
615
616 int
atom_to_unique_offset(int iatom) const617 Molecule::atom_to_unique_offset(int iatom) const
618 {
619 int iuniq = atom_to_uniq_[iatom];
620 int nequiv = nequiv_[iuniq];
621 for (int i=0; i<nequiv; i++) {
622 if (equiv_[iuniq][i] == iatom) return i;
623 }
624 ExEnv::errn() << "Molecule::atom_to_unique_offset: internal error"
625 << endl;
626 return -1;
627 }
628
629 void
set_point_group(const Ref<PointGroup> & ppg,double tol)630 Molecule::set_point_group(const Ref<PointGroup>&ppg, double tol)
631 {
632 ExEnv::out0() << indent
633 << "Molecule: setting point group to " << ppg->symbol()
634 << endl;
635 pg_ = new PointGroup(*ppg.pointer());
636
637 double r[3];
638 for (int i=0; i<3; i++) {
639 r[i] = -pg_->origin()[i];
640 pg_->origin()[i] = 0;
641 }
642 translate(r);
643
644 clear_symmetry_info();
645 init_symmetry_info();
646
647 cleanup_molecule(tol);
648 }
649
650 Ref<PointGroup>
point_group() const651 Molecule::point_group() const
652 {
653 return pg_;
654 }
655
656 SCVector3
center_of_mass() const657 Molecule::center_of_mass() const
658 {
659 SCVector3 ret;
660 double M;
661
662 ret = 0.0;
663 M = 0.0;
664
665 for (int i=0; i < natom(); i++) {
666 double m = mass(i);
667 ret += m * SCVector3(r(i));
668 M += m;
669 }
670
671 ret *= 1.0/M;
672
673 return ret;
674 }
675
676 SCVector3
center_of_charge() const677 Molecule::center_of_charge() const
678 {
679 SCVector3 ret;
680 double C;
681
682 ret = 0.0;
683 C = 0.0;
684
685 for (int i=0; i < natom(); i++) {
686 double c = charge(i);
687 ret += c * SCVector3(r(i));
688 C += c;
689 }
690
691 ret *= 1.0/C;
692
693 return ret;
694 }
695
696
697 double
nuclear_repulsion_energy()698 Molecule::nuclear_repulsion_energy()
699 {
700 double e=0.0;
701
702 // non_q non_q terms
703 for (int ii=1; ii < non_q_atoms_.size(); ii++) {
704 int i = non_q_atoms_[ii];
705 SCVector3 ai(r(i));
706 double Zi = charge(i);
707
708 for (int jj=0; jj < ii; jj++) {
709 int j = non_q_atoms_[jj];
710 SCVector3 aj(r(j));
711 e += Zi * charge(j) / ai.dist(aj);
712 }
713 }
714
715 // non_q q terms
716 for (int ii=0; ii < q_atoms_.size(); ii++) {
717 int i = q_atoms_[ii];
718 SCVector3 ai(r(i));
719 double Zi = charge(i);
720
721 for (int jj=0; jj < non_q_atoms_.size(); jj++) {
722 int j = non_q_atoms_[jj];
723 SCVector3 aj(r(j));
724 e += Zi * charge(j) / ai.dist(aj);
725 }
726 }
727
728 if (include_qq_) {
729 // q q terms
730 for (int ii=1; ii < q_atoms_.size(); ii++) {
731 int i = q_atoms_[ii];
732 SCVector3 ai(r(i));
733 double Zi = charge(i);
734
735 for (int jj=0; jj < ii; jj++) {
736 int j = q_atoms_[jj];
737 SCVector3 aj(r(j));
738 e += Zi * charge(j) / ai.dist(aj);
739 }
740 }
741 }
742
743 return e;
744 }
745
746 void
nuclear_repulsion_1der(int center,double xyz[3])747 Molecule::nuclear_repulsion_1der(int center, double xyz[3])
748 {
749 int i,j,k;
750 double rd[3],r2;
751 double factor;
752
753 xyz[0] = 0.0;
754 xyz[1] = 0.0;
755 xyz[2] = 0.0;
756
757 SCVector3 r_center(r(center));
758 double Z_center = charge(center);
759 bool center_is_Q = (atom_symbol(center) == "Q");
760
761 // this handles center = Q or non_Q and atom = non_Q
762 for (int ii=0; ii < non_q_atoms_.size(); ii++) {
763 int i = non_q_atoms_[ii];
764 if (i == center) continue;
765 SCVector3 r_i(r(i));
766
767 r2 = 0.0;
768 for (k=0; k < 3; k++) {
769 rd[k] = r_center[k] - r_i[k];
770 r2 += rd[k]*rd[k];
771 }
772 factor = - Z_center * charge(i) * pow(r2,-1.5);
773 for (k=0; k<3; k++) {
774 xyz[k] += factor * rd[k];
775 }
776 }
777
778 // this handles center = Q or non_Q and atom = Q
779 for (int ii=0; ii < q_atoms_.size(); ii++) {
780 int i = q_atoms_[ii];
781 if (i == center || (!include_qq_ && center_is_Q)) continue;
782 SCVector3 r_i(r(i));
783
784 r2 = 0.0;
785 for (k=0; k < 3; k++) {
786 rd[k] = r_center[k] - r_i[k];
787 r2 += rd[k]*rd[k];
788 }
789 factor = - Z_center * charge(i) * pow(r2,-1.5);
790 for (k=0; k<3; k++) {
791 xyz[k] += factor * rd[k];
792 }
793 }
794 }
795
796 void
nuclear_charge_efield(const double * charges,const double * position,double * efield)797 Molecule::nuclear_charge_efield(const double *charges,
798 const double *position, double *efield)
799 {
800 double tmp;
801 double rd[3];
802
803 for (int i=0; i<3; i++) efield[i] = 0.0;
804
805 if (include_q_) {
806 for (int i=0; i<natoms_; i++) {
807 SCVector3 a(r(i));
808 tmp = 0.0;
809 for (int j=0; j<3; j++) {
810 rd[j] = position[j] - a[j];
811 tmp += rd[j]*rd[j];
812 }
813 tmp = charges[i]/(tmp*sqrt(tmp));
814 for (int j=0; j<3; j++) {
815 efield[j] += rd[j] * tmp;
816 }
817 }
818 }
819 else {
820 for (int ii=0; ii<non_q_atoms_.size(); ii++) {
821 int i = non_q_atoms_[ii];
822 SCVector3 a(r(i));
823 tmp = 0.0;
824 for (int j=0; j<3; j++) {
825 rd[j] = position[j] - a[j];
826 tmp += rd[j]*rd[j];
827 }
828 tmp = charges[i]/(tmp*sqrt(tmp));
829 for (int j=0; j<3; j++) {
830 efield[j] += rd[j] * tmp;
831 }
832 }
833 }
834 }
835
836 void
nuclear_efield(const double * position,double * efield)837 Molecule::nuclear_efield(const double *position, double *efield)
838 {
839 double tmp;
840 double rd[3];
841
842 for (int i=0; i<3; i++) efield[i] = 0.0;
843
844 if (include_q_) {
845 for (int i=0; i<natoms_; i++) {
846 SCVector3 a(r(i));
847 tmp = 0.0;
848 for (int j=0; j<3; j++) {
849 rd[j] = position[j] - a[j];
850 tmp += rd[j]*rd[j];
851 }
852 tmp = charge(i)/(tmp*sqrt(tmp));
853 for (int j=0; j<3; j++) {
854 efield[j] += rd[j] * tmp;
855 }
856 }
857 }
858 else {
859 for (int ii=0; ii<non_q_atoms_.size(); ii++) {
860 int i = non_q_atoms_[ii];
861 SCVector3 a(r(i));
862 tmp = 0.0;
863 for (int j=0; j<3; j++) {
864 rd[j] = position[j] - a[j];
865 tmp += rd[j]*rd[j];
866 }
867 tmp = charge(i)/(tmp*sqrt(tmp));
868 for (int j=0; j<3; j++) {
869 efield[j] += rd[j] * tmp;
870 }
871 }
872 }
873 }
874
875 int
atom_at_position(double * v,double tol) const876 Molecule::atom_at_position(double *v, double tol) const
877 {
878 SCVector3 p(v);
879 for (int i=0; i < natom(); i++) {
880 SCVector3 ai(r(i));
881 if (p.dist(ai) < tol) return i;
882 }
883 return -1;
884 }
885
886 void
symmetrize(const Ref<PointGroup> & pg,double tol)887 Molecule::symmetrize(const Ref<PointGroup> &pg, double tol)
888 {
889 pg_ = new PointGroup(pg);
890
891 // translate to the origin of the symmetry frame
892 double r[3];
893 for (int i=0; i<3; i++) {
894 r[i] = -pg_->origin()[i];
895 pg_->origin()[i] = 0;
896 }
897 translate(r);
898
899 symmetrize(tol);
900 }
901
902 // We are given a molecule which may or may not have just the symmetry
903 // distinct atoms in it. We have to go through the existing set of atoms,
904 // perform each symmetry operation in the point group on each of them, and
905 // then add the new atom if it isn't in the list already
906
907 void
symmetrize(double tol)908 Molecule::symmetrize(double tol)
909 {
910 // if molecule is c1, don't do anything
911 if (!strcmp(this->point_group()->symbol(),"c1")) {
912 init_symmetry_info();
913 return;
914 }
915
916 clear_symmetry_info();
917
918 Molecule *newmol = new Molecule(*this);
919
920 CharacterTable ct = this->point_group()->char_table();
921
922 SCVector3 np;
923 SymmetryOperation so;
924
925 for (int i=0; i < natom(); i++) {
926 SCVector3 ac(r(i));
927
928 for (int g=0; g < ct.order(); g++) {
929 so = ct.symm_operation(g);
930 for (int ii=0; ii < 3; ii++) {
931 np[ii]=0;
932 for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * ac[jj];
933 }
934
935 int atom = newmol->atom_at_position(np.data(), tol);
936 if (atom < 0) {
937 newmol->add_atom(Z_[i],np[0],np[1],np[2],label(i));
938 }
939 else {
940 if (Z(i) != newmol->Z(atom)
941 || fabs(mass(i)-newmol->mass(atom))>1.0e-10) {
942 throw ToleranceExceeded("symmetrize: atom mismatch",
943 __FILE__, __LINE__,
944 1.0e-10, fabs(mass(i)-newmol->mass(atom)),
945 class_desc());
946 }
947 }
948 }
949 }
950
951 Ref<Units> saved_units = geometry_units_;
952 *this = *newmol;
953 geometry_units_ = saved_units;
954 delete newmol;
955
956 init_symmetry_info();
957 }
958
959 void
translate(const double * r)960 Molecule::translate(const double *r)
961 {
962 for (int i=0; i < natom(); i++) {
963 r_[i][0] += r[0];
964 r_[i][1] += r[1];
965 r_[i][2] += r[2];
966 }
967 }
968
969 // move the molecule to the center of mass
970 void
move_to_com()971 Molecule::move_to_com()
972 {
973 SCVector3 com = -center_of_mass();
974 translate(com.data());
975 }
976
977 // find the 3 principal coordinate axes, and rotate the molecule to be
978 // aligned along them. also rotate the symmetry frame contained in point_group
979 void
transform_to_principal_axes(int trans_frame)980 Molecule::transform_to_principal_axes(int trans_frame)
981 {
982 // mol_move_to_com(mol);
983
984 double *inert[3], inert_dat[9], *evecs[3], evecs_dat[9];
985 double evals[3];
986
987 int i,j,k;
988 for (i=0; i < 3; i++) {
989 inert[i] = &inert_dat[i*3];
990 evecs[i] = &evecs_dat[i*3];
991 }
992 memset(inert_dat,'\0',sizeof(double)*9);
993 memset(evecs_dat,'\0',sizeof(double)*9);
994
995 for (i=0; i < natom(); i++) {
996 SCVector3 ac(r(i));
997 double m=mass(i);
998 inert[0][0] += m * (ac[1]*ac[1] + ac[2]*ac[2]);
999 inert[1][0] -= m * ac[0]*ac[1];
1000 inert[1][1] += m * (ac[0]*ac[0] + ac[2]*ac[2]);
1001 inert[2][0] -= m * ac[0]*ac[2];
1002 inert[2][1] -= m * ac[1]*ac[2];
1003 inert[2][2] += m * (ac[0]*ac[0] + ac[1]*ac[1]);
1004 }
1005
1006 inert[0][1] = inert[1][0];
1007 inert[0][2] = inert[2][0];
1008 inert[1][2] = inert[2][1];
1009
1010 // cleanup inert
1011 for (i=0; i < 3; i++) {
1012 for (int j=0; j <= i; j++) {
1013 if (fabs(inert[i][j]) < 1.0e-5) {
1014 inert[i][j]=inert[j][i]=0.0;
1015 }
1016 }
1017 }
1018
1019 cmat_diag(inert, evals, evecs, 3, 1, 1e-14);
1020
1021 // cleanup evecs
1022 for (i=0; i < 3; i++) {
1023 for (int j=0; j < 3; j++) {
1024 if (fabs(evecs[i][j]) < 1.0e-5) {
1025 evecs[i][j]=0.0;
1026 }
1027 }
1028 }
1029
1030 for (i=0; i<natom(); i++) {
1031 double a[3];
1032 a[0] = r(i,0); a[1] = r(i,1); a[2] = r(i,2);
1033 for (j=0; j<3; j++) {
1034 double e = 0.0;
1035 for (k=0; k<3; k++) {
1036 e += a[k] * evecs[k][j];
1037 }
1038 r_[i][j] = e;
1039 }
1040 }
1041
1042 if (!trans_frame) return;
1043
1044 SymmetryOperation tso=point_group()->symm_frame();
1045
1046 for (i=0; i < 3; i++) {
1047 for (int j=0; j < 3; j++) {
1048 double t=0;
1049 for (int k=0; k < 3; k++) t += tso[k][j]*evecs[k][i];
1050 pg_->symm_frame()[i][j] = t;
1051 }
1052 }
1053 }
1054
1055 void
transform_to_symmetry_frame()1056 Molecule::transform_to_symmetry_frame()
1057 {
1058 int i,j,k;
1059 double t[3][3];
1060
1061 SymmetryOperation tso=point_group()->symm_frame();
1062
1063 for (i=0; i<3; i++) {
1064 for (j=0; j<3; j++) {
1065 t[i][j] = tso[i][j];
1066 }
1067 }
1068
1069 for (i=0; i<natom(); i++) {
1070 double a[3];
1071 a[0] = r(i,0); a[1] = r(i,1); a[2] = r(i,2);
1072 for (j=0; j<3; j++) {
1073 double e = 0.0;
1074 for (k=0; k<3; k++) {
1075 e += a[k] * t[k][j];
1076 }
1077 r_[i][j] = e;
1078 }
1079 }
1080
1081 for (i=0; i<3; i++) {
1082 for (j=0; j<3; j++) {
1083 double e=0;
1084 for (k=0; k<3; k++) e += tso[k][j]*t[k][i];
1085 pg_->symm_frame()[i][j] = e;
1086 }
1087 }
1088 }
1089
1090 // given a molecule, make sure that equivalent centers have coordinates
1091 // that really map into each other
1092
1093 void
cleanup_molecule(double tol)1094 Molecule::cleanup_molecule(double tol)
1095 {
1096 // if symmetry is c1, do nothing else
1097 if (!strcmp(point_group()->symbol(),"c1")) return;
1098
1099 int i;
1100 SCVector3 up,np,ap;
1101 SymmetryOperation so;
1102 CharacterTable ct = point_group()->char_table();
1103
1104 // first clean up the unique atoms by replacing each coordinate with the
1105 // average of coordinates obtained by applying all symmetry operations to
1106 // the original atom, iff the new atom ends up near the original atom
1107 for (i=0; i < nunique(); i++) {
1108 // up will store the original coordinates of unique atom i
1109 up = r(unique(i));
1110 // ap will hold the average coordinate (times the number of coordinates)
1111 // initialize it to the E result
1112 ap = up;
1113 int ncoor = 1;
1114 // loop through all sym ops except E
1115 for (int g=1; g < ct.order(); g++) {
1116 so = ct.symm_operation(g);
1117 for (int ii=0; ii < 3; ii++) {
1118 np[ii]=0;
1119 for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * up[jj];
1120 }
1121 if (np.dist(up) < 0.1) {
1122 for (int jj=0; jj < 3; jj++) ap[jj] += np[jj];
1123 ncoor++;
1124 }
1125 }
1126 // replace the unique coordinate with the average coordinate
1127 r_[unique(i)][0] = ap[0] / ncoor;
1128 r_[unique(i)][1] = ap[1] / ncoor;
1129 r_[unique(i)][2] = ap[2] / ncoor;
1130 }
1131
1132 // find the atoms equivalent to each unique atom and eliminate
1133 // numerical errors that may be in the equivalent atom's coordinates
1134
1135 // loop through unique atoms
1136 for (i=0; i < nunique(); i++) {
1137 // up will store the coordinates of unique atom i
1138 up = r(unique(i));
1139
1140 // loop through all sym ops except E
1141 for (int g=1; g < ct.order(); g++) {
1142 so = ct.symm_operation(g);
1143 for (int ii=0; ii < 3; ii++) {
1144 np[ii]=0;
1145 for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * up[jj];
1146 }
1147
1148 // loop through equivalent atoms
1149 int found = 0;
1150 for (int j=0; j < natom(); j++) {
1151 // see if j is generated from i
1152 if (np.dist(SCVector3(r(j))) < tol) {
1153 r_[j][0] = np[0];
1154 r_[j][1] = np[1];
1155 r_[j][2] = np[2];
1156 found = 1;
1157 }
1158 }
1159 if (!found) {
1160 SCException ex("cleanup: couldn't find atom",
1161 __FILE__, __LINE__, class_desc());
1162 try {
1163 ex.elaborate()
1164 << "couldn't find atom at " << np << endl
1165 << "transforming uniq atom " << i << " at " << up << endl
1166 << "with symmetry op " << g << ":" << endl;
1167 so.print(ex.elaborate());
1168 }
1169 catch (...) {}
1170 throw ex;
1171 }
1172 }
1173 }
1174
1175 }
1176
1177 ///////////////////////////////////////////////////////////////////
1178 // Compute the principal axes and the principal moments of inertia
1179 ///////////////////////////////////////////////////////////////////
1180
1181 void
principal_moments_of_inertia(double * evals,double ** evecs) const1182 Molecule::principal_moments_of_inertia(double *evals, double **evecs) const
1183 {
1184
1185 // The principal moments of inertia are computed in amu*angstrom^2
1186 // evals: principal moments of inertia
1187 // evecs: principal axes (optional argument)
1188
1189 Ref<Units> units = new Units("angstroms * angstroms");
1190 double au_to_angs = units->from_atomic_units();
1191
1192 double *inert[3]; // inertia tensor
1193
1194 int i, j;
1195 int delete_evecs = 0;
1196
1197 // (allocate and) initialize evecs, evals, and inert
1198 if (!evecs) {
1199 evecs = new double*[3];
1200 for (i=0; i<3; i++) evecs[i] = new double[3];
1201 delete_evecs = 1;
1202 }
1203 for (i=0; i<3; i++) {
1204 inert[i] = new double[3];
1205 memset(inert[i],'\0',sizeof(double)*3);
1206 memset(evecs[i],'\0',sizeof(double)*3);
1207 }
1208 memset(evals,'\0',sizeof(double)*3);
1209
1210 SCVector3 com = center_of_mass();
1211
1212 // compute inertia tensor
1213 SCVector3 ac;
1214 for (i=0; i<natom(); i++) {
1215 ac = r(i);
1216 // compute moments of inertia wrt center of mass
1217 for (j=0; j<3; j++) ac(j) -= com(j);
1218 double m=au_to_angs*mass(i);
1219 inert[0][0] += m * (ac[1]*ac[1] + ac[2]*ac[2]);
1220 inert[1][0] -= m * ac[0]*ac[1];
1221 inert[1][1] += m * (ac[0]*ac[0] + ac[2]*ac[2]);
1222 inert[2][0] -= m * ac[0]*ac[2];
1223 inert[2][1] -= m * ac[1]*ac[2];
1224 inert[2][2] += m * (ac[0]*ac[0] + ac[1]*ac[1]);
1225 }
1226 inert[0][1] = inert[1][0];
1227 inert[0][2] = inert[2][0];
1228 inert[1][2] = inert[2][1];
1229
1230 cmat_diag(inert, evals, evecs, 3, 1, 1e-14);
1231
1232 if (delete_evecs) {
1233 for (i=0; i<3; i++) delete[] evecs[i];
1234 delete[] evecs;
1235 }
1236 for (i=0; i<3; i++) {
1237 delete[] inert[i];
1238 }
1239 }
1240
1241 int
n_core_electrons()1242 Molecule::n_core_electrons()
1243 {
1244 int i,n=0;
1245 for (i=0; i<natom(); i++) {
1246 if (charge(i) == 0.0) continue;
1247 int z = Z_[i];
1248 if (z > 2) n += 2;
1249 if (z > 10) n += 8;
1250 if (z > 18) n += 8;
1251 if (z > 30) n += 10;
1252 if (z > 36) n += 8;
1253 if (z > 48) n += 10;
1254 if (z > 54) n += 8;
1255 if (z > 72) {
1256 throw LimitExceeded<int>("n_core_electrons: atomic number too large",
1257 __FILE__, __LINE__, 72, z, class_desc());
1258 }
1259 }
1260 return n;
1261 }
1262
1263 int
max_z()1264 Molecule::max_z()
1265 {
1266 int i, maxz=0;
1267 for (i=0; i<natom(); i++) {
1268 int z = Z_[i];
1269 if (z>maxz) maxz = z;
1270 }
1271 return maxz;
1272 }
1273
1274 void
read_pdb(const char * filename)1275 Molecule::read_pdb(const char *filename)
1276 {
1277 clear();
1278 ifstream in(filename);
1279 Ref<Units> units = new Units("angstrom");
1280 while (in.good()) {
1281 const int max_line = 80;
1282 char line[max_line];
1283 in.getline(line,max_line);
1284 char *endofline = (char*) memchr(line, 0, max_line);
1285 if (endofline) memset(endofline, ' ', &line[max_line-1] - endofline);
1286 if (!in.good()) break;
1287 if (strncmp(line,"ATOM ",6) == 0
1288 ||strncmp(line,"HETATM",6) == 0) {
1289
1290 char element[3];
1291 strncpy(element,&line[76],2); element[2] = '\0';
1292 char name[5];
1293 strncpy(name,&line[12],4); name[4] = '\0';
1294 if (element[0]==' '&&element[1]==' ') {
1295 // no element was given so get the element from the atom name
1296 if (name[0]!=' '&&name[3]!=' ') {
1297 // some of the atom label may have been
1298 // pushed over into the element fields
1299 // so check the residue
1300 char resName[4];
1301 strncpy(resName,&line[17],3); resName[3] = '\0';
1302 if (strncmp(line,"ATOM ",6)==0&&(0
1303 ||strcmp(resName,"ALA")==0||strcmp(resName,"A ")==0
1304 ||strcmp(resName,"ARG")==0||strcmp(resName,"R ")==0
1305 ||strcmp(resName,"ASN")==0||strcmp(resName,"N ")==0
1306 ||strcmp(resName,"ASP")==0||strcmp(resName,"D ")==0
1307 ||strcmp(resName,"ASX")==0||strcmp(resName,"B ")==0
1308 ||strcmp(resName,"CYS")==0||strcmp(resName,"C ")==0
1309 ||strcmp(resName,"GLN")==0||strcmp(resName,"Q ")==0
1310 ||strcmp(resName,"GLU")==0||strcmp(resName,"E ")==0
1311 ||strcmp(resName,"GLX")==0||strcmp(resName,"Z ")==0
1312 ||strcmp(resName,"GLY")==0||strcmp(resName,"G ")==0
1313 ||strcmp(resName,"HIS")==0||strcmp(resName,"H ")==0
1314 ||strcmp(resName,"ILE")==0||strcmp(resName,"I ")==0
1315 ||strcmp(resName,"LEU")==0||strcmp(resName,"L ")==0
1316 ||strcmp(resName,"LYS")==0||strcmp(resName,"K ")==0
1317 ||strcmp(resName,"MET")==0||strcmp(resName,"M ")==0
1318 ||strcmp(resName,"PHE")==0||strcmp(resName,"F ")==0
1319 ||strcmp(resName,"PRO")==0||strcmp(resName,"P ")==0
1320 ||strcmp(resName,"SER")==0||strcmp(resName,"S ")==0
1321 ||strcmp(resName,"THR")==0||strcmp(resName,"T ")==0
1322 ||strcmp(resName,"TRP")==0||strcmp(resName,"W ")==0
1323 ||strcmp(resName,"TYR")==0||strcmp(resName,"Y ")==0
1324 ||strcmp(resName,"VAL")==0||strcmp(resName,"V ")==0
1325 ||strcmp(resName,"A ")==0
1326 ||strcmp(resName,"+A ")==0
1327 ||strcmp(resName,"C ")==0
1328 ||strcmp(resName,"+C ")==0
1329 ||strcmp(resName,"G ")==0
1330 ||strcmp(resName,"+G ")==0
1331 ||strcmp(resName,"I ")==0
1332 ||strcmp(resName,"+I ")==0
1333 ||strcmp(resName,"T ")==0
1334 ||strcmp(resName,"+T ")==0
1335 ||strcmp(resName,"U ")==0
1336 ||strcmp(resName,"+U ")==0
1337 )) {
1338 // there no two letter elements for these cases
1339 element[0] = name[0];
1340 element[1] = '\0';
1341 }
1342 }
1343 else {
1344 strncpy(element,name,2); element[2] = '\0';
1345 }
1346 }
1347 if (element[0] == ' ') {
1348 element[0] = element[1];
1349 element[1] = '\0';
1350 }
1351
1352 int Z = atominfo_->string_to_Z(element);
1353
1354 char field[9];
1355 strncpy(field,&line[30],8); field[8] = '\0';
1356 double x = atof(field);
1357 strncpy(field,&line[38],8); field[8] = '\0';
1358 double y = atof(field);
1359 strncpy(field,&line[46],8); field[8] = '\0';
1360 double z = atof(field);
1361 add_atom(Z,
1362 x*units->to_atomic_units(),
1363 y*units->to_atomic_units(),
1364 z*units->to_atomic_units());
1365 }
1366 else {
1367 // skip to next record
1368 }
1369 }
1370 }
1371
1372 void
print_pdb(ostream & os,char * title) const1373 Molecule::print_pdb(ostream& os, char *title) const
1374 {
1375 Ref<Units> u = new Units("angstrom");
1376 double bohr = u->from_atomic_units();
1377
1378 if (title)
1379 os << scprintf("%-10s%-60s\n","COMPND",title);
1380 else
1381 os << scprintf("%-10s%-60s\n","COMPND","Title");
1382
1383 if (title)
1384 os << scprintf("REMARK %s\n", title);
1385
1386 int i;
1387 for (i=0; i < natom(); i++) {
1388 char symb[4];
1389 std::string symbol(atom_symbol(i));
1390 sprintf(symb,"%s1",symbol.c_str());
1391
1392 os << scprintf(
1393 "HETATM%5d %-3s UNK %5d %8.3f%8.3f%8.3f 0.00 0.00 0\n",
1394 i+1, symb, 0, r(i,0)*bohr, r(i,1)*bohr, r(i,2)*bohr);
1395 }
1396
1397 for (i=0; i < natom(); i++) {
1398 double at_rad_i = atominfo_->atomic_radius(Z_[i]);
1399 SCVector3 ai(r(i));
1400
1401 os << scprintf("CONECT%5d",i+1);
1402
1403 for (int j=0; j < natom(); j++) {
1404
1405 if (j==i) continue;
1406
1407 double at_rad_j = atominfo_->atomic_radius(Z_[j]);
1408 SCVector3 aj(r(j));
1409
1410 if (ai.dist(aj) < 1.1*(at_rad_i+at_rad_j))
1411 os << scprintf("%5d",j+1);
1412 }
1413
1414 os << endl;
1415 }
1416
1417 os << "END" << endl;
1418 os.flush();
1419 }
1420
1421 double
mass(int atom) const1422 Molecule::mass(int atom) const
1423 {
1424 if (!mass_ || mass_[atom] == 0) {
1425 return atominfo_->mass(Z_[atom]);
1426 }
1427 return mass_[atom];
1428 }
1429
1430 const char *
label(int atom) const1431 Molecule::label(int atom) const
1432 {
1433 if (!labels_) return 0;
1434 return labels_[atom];
1435 }
1436
1437 std::string
atom_name(int iatom) const1438 Molecule::atom_name(int iatom) const
1439 {
1440 return atominfo_->name(Z_[iatom]);
1441 }
1442
1443 std::string
atom_symbol(int iatom) const1444 Molecule::atom_symbol(int iatom) const
1445 {
1446 return atominfo_->symbol(Z_[iatom]);
1447 }
1448
1449 /////////////////////////////////////////////////////////////////////////////
1450
1451 // Local Variables:
1452 // mode: c++
1453 // c-file-style: "CLJ"
1454 // End:
1455