1 // PROJECT.CPP
2
3 // Copyright (C) 1998 Tommi Hassinen.
4
5 // This package is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation; either version 2 of the License, or
8 // (at your option) any later version.
9
10 // This package is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14
15 // You should have received a copy of the GNU General Public License
16 // along with this package; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 /*################################################################################################*/
20
21 #include "project.h" // config.h is here -> we get ENABLE-macros here...
22
23 #include <ghemical/libghemicaldefine.h>
24
25 #include <ghemical/v3d.h>
26 #include <ghemical/atom.h>
27 #include <ghemical/bond.h>
28
29 #include <ghemical/seqbuild.h>
30
31 #include <ghemical/engine.h>
32
33 #include <ghemical/eng1_qm.h>
34 #include <ghemical/eng1_mm.h>
35 #include <ghemical/eng1_sf.h>
36
37 #include <ghemical/eng2_qm_mm.h>
38
39 #include <ghemical/geomopt.h>
40 #include <ghemical/intcrd.h>
41
42 #include <ghemical/pop_ana.h>
43
44 #include <ghemical/search.h>
45 #include <ghemical/utility.h>
46
47 #include <ghemical/notice.h>
48
49 #include "appdefine.h"
50
51 #include "custom_app.h"
52 #include "custom_camera.h"
53 #include "custom_lights.h"
54
55 #include "ogl_plane.h"
56 #include "ogl_surface.h"
57
58 #include "color.h"
59
60 #include "filetrans.h"
61
62 #include "local_i18n.h"
63
64 #include <cstring>
65 #include <fstream>
66 #include <sstream>
67 #include <iomanip>
68 #include <algorithm>
69 using namespace std;
70
71 /*################################################################################################*/
72
custom_transformer_client(void)73 custom_transformer_client::custom_transformer_client(void) :
74 ogl_transformer_client()
75 {
76 tc_object_ref = NULL;
77 tc_local_object = false;
78 }
79
~custom_transformer_client(void)80 custom_transformer_client::~custom_transformer_client(void)
81 {
82 }
83
84 /*################################################################################################*/
85
86 ogl_dummy_object * project::selected_object = NULL;
87
88 const char project::appversion[16] = APPVERSION;
89 char project::appdata_path[256] = "appdata_path_is_not_yet_set_by_main_program";
90
91 bool project::background_job_running = false;
92
93 iGLu project::list_counter = 1; // zero is not a valid display list number...
94
95 color_mode_element project::cm_element = color_mode_element();
96 color_mode_secstruct project::cm_secstruct = color_mode_secstruct();
97 color_mode_hydphob project::cm_hydphob = color_mode_hydphob();
98
project(void)99 project::project(void) :
100 custom_transformer_client(),
101 model()
102 {
103 project_path = NULL;
104 project_filename = NULL;
105
106 SetDefaultProjectFileName();
107
108 selected_object = NULL; // always re-initialize this!!!
109
110 camera_counter = 1;
111 object_counter = 1;
112
113 mt_a1 = mt_a2 = mt_a3 = NULL;
114
115 importpdb_mdata = NULL; // temporary?!?!?!
116 }
117
~project(void)118 project::~project(void)
119 {
120 selected_object = NULL;
121
122 if (object_vector.size() != 0) assertion_failed(__FILE__, __LINE__, "object_vector.size() != 0");
123
124 if (plotting_view_vector.size() != 0) assertion_failed(__FILE__, __LINE__, "plotting_view_vector.size() != 0");
125
126 if (graphics_view_vector.size() != 0) assertion_failed(__FILE__, __LINE__, "graphics_view_vector.size() != 0");
127
128 if (bond_list.size() != 0) assertion_failed(__FILE__, __LINE__, "bond_list.size() != 0");
129
130 if (atom_list.size() != 0) assertion_failed(__FILE__, __LINE__, "atom_list.size() != 0");
131
132 if (project_path != NULL)
133 {
134 delete[] project_path;
135 project_path = NULL;
136 }
137
138 if (project_filename != NULL)
139 {
140 delete[] project_filename;
141 project_filename = NULL;
142 }
143 }
144
ClearAll(void)145 void project::ClearAll(void)
146 {
147 selected_object = NULL;
148
149 while (object_vector.size() != 0)
150 {
151 ogl_smart_object * ref;
152 ref = object_vector.back();
153 RemoveObject(ref);
154 }
155
156 while (plotting_view_vector.size() != 0)
157 {
158 base_wcl * ref;
159 ref = plotting_view_vector.back();
160 RemovePlottingClient(ref);
161 }
162
163 while (graphics_view_vector.size() != 0)
164 {
165 oglview_wcl * ref;
166 ref = graphics_view_vector.back();
167 RemoveGraphicsClient(ref, true);
168 }
169
170 while (bond_list.size() != 0)
171 {
172 iter_bl itB = GetBondsBegin();
173 RemoveBond(itB);
174 }
175
176 while (atom_list.size() != 0)
177 {
178 iter_al itA = GetAtomsBegin();
179 RemoveAtom(itA);
180 }
181 }
182
GetProjectFileNameExtension(void)183 const char * project::GetProjectFileNameExtension(void)
184 {
185 static const char ext[] = "gpr";
186 return ext;
187 }
188
SetProjectPath(const char * path)189 void project::SetProjectPath(const char * path)
190 {
191 if (project_path != NULL) delete[] project_path;
192
193 project_path = new char[strlen(path) + 1];
194 strcpy(project_path, path);
195 }
196
SetProjectFileName(const char * filename)197 void project::SetProjectFileName(const char * filename)
198 {
199 if (project_filename != NULL) delete[] project_filename;
200
201 project_filename = new char[strlen(filename) + 1];
202 strcpy(project_filename, filename);
203 }
204
SetDefaultProjectFileName(void)205 void project::SetDefaultProjectFileName(void)
206 {
207 static i32s id_counter = 1;
208
209 ostringstream str;
210 str << _("untitled") << setw(2) << setfill('0') << id_counter++ << ends;
211
212 SetProjectFileName(str.str().c_str());
213 }
214
ParseProjectFileNameAndPath(const char * string)215 void project::ParseProjectFileNameAndPath(const char * string)
216 {
217 char * localstring1 = new char[strlen(string) + 1];
218 strcpy(localstring1, string);
219
220 i32s lastdir = NOT_DEFINED;
221 for (i32s n1 = 0;n1 < (i32s) strlen(localstring1);n1++)
222 {
223 if (localstring1[n1] == DIR_SEPARATOR) lastdir = n1;
224 }
225
226 char * localstring2 = localstring1;
227
228 // set project_path if needed...
229 // set project_path if needed...
230 // set project_path if needed...
231
232 if (lastdir != NOT_DEFINED)
233 {
234 localstring2[lastdir] = 0;
235 SetProjectPath(localstring2);
236
237 localstring2 = & localstring2[lastdir + 1];
238 }
239
240 // and now set project_filename, without extension...
241 // and now set project_filename, without extension...
242 // and now set project_filename, without extension...
243
244 i32s lastext = NOT_DEFINED;
245 for (i32s n1 = 0;n1 < (i32s) strlen(localstring2);n1++)
246 {
247 if (localstring2[n1] == EXT_SEPARATOR) lastext = n1;
248 }
249
250 if (lastext != NOT_DEFINED)
251 {
252 // this only removes an extension if it matches *our* extension,
253 // which makes problems for imported files e.g. nh3.mol.mmg1p (!) instead of nh3.mol or nh3.mm1gp
254 // char * localstring3 = & localstring2[lastext + 1];
255 // bool has_extension = !strcmp(localstring3, GetProjectFileNameExtension());
256 // if (has_extension) localstring2[lastext] = 0;
257
258 // use this instead:
259 // ^^^^^^^^^^^^^^^^^
260
261 localstring2[lastext] = 0;
262 }
263
264 SetProjectFileName(localstring2);
265
266 delete[] localstring1;
267 }
268
GetProjectFileName(char * buffer,int buffer_size,bool with_extension)269 void project::GetProjectFileName(char * buffer, int buffer_size, bool with_extension)
270 {
271 ostringstream ostr;
272
273 ostr << project_filename;
274 if (with_extension) ostr << EXT_SEPARATOR << GetProjectFileNameExtension();
275 ostr << ends;
276
277 if (strlen(ostr.str().c_str()) + 1 >= buffer_size)
278 {
279 assertion_failed(__FILE__, __LINE__, "buffer overflow!");
280 }
281
282 strcpy(buffer, ostr.str().c_str());
283 }
284
GetFullProjectFileName(char * buffer,int buffer_size)285 void project::GetFullProjectFileName(char * buffer, int buffer_size)
286 {
287 ostringstream ostr;
288
289 if (project_path != NULL) ostr << project_path << DIR_SEPARATOR;
290 ostr << project_filename << EXT_SEPARATOR << GetProjectFileNameExtension() << ends;
291
292 if (strlen(ostr.str().c_str()) + 1 >= buffer_size)
293 {
294 assertion_failed(__FILE__, __LINE__, "buffer overflow!");
295 }
296
297 strcpy(buffer, ostr.str().c_str());
298 }
299
300 /*##############################################*/
301 /*##############################################*/
302
303 #ifdef ENABLE_OPENBABEL
304
ImportFile(const char * filename,int index)305 bool project::ImportFile(const char * filename, int index)
306 {
307 ifstream ifile;
308 ostringstream intermed;
309 file_trans translator;
310
311 // store the current (numeric) locale into my_num_locale,
312 // and switch into the "C" numeric locale...
313
314 static char my_num_locale[32] = "C";
315 strcpy(my_num_locale, setlocale(LC_NUMERIC, NULL));
316 setlocale(LC_NUMERIC, "C");
317
318 if (index == 0) // Automatic detection
319 {
320 if (!translator.CanImport(filename))
321 {
322 ErrorMessage(_("Cannot import that file type."));
323 return false;
324 }
325
326 ifile.open(filename, ios::in);
327 translator.Import(filename, ifile, intermed);
328 ifile.close();
329 }
330 else // By type picked by the user
331 {
332 ifile.open(filename, ios::in);
333 translator.Import(filename, index - 1, ifile, intermed);
334 ifile.close();
335 }
336
337 istringstream interInput(intermed.str());
338 bool retval = ReadGPR((* this), interInput, false);
339
340 // change the original locale back...
341 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
342 setlocale(LC_NUMERIC, my_num_locale);
343
344 return retval;
345 }
346
ExportFile(const char * filename,int index)347 bool project::ExportFile(const char * filename, int index)
348 {
349 ofstream ofile;
350 stringstream intermed;
351 file_trans translator;
352
353 // store the current (numeric) locale into my_num_locale,
354 // and switch into the "C" numeric locale...
355
356 static char my_num_locale[32] = "C";
357 strcpy(my_num_locale, setlocale(LC_NUMERIC, NULL));
358 setlocale(LC_NUMERIC, "C");
359
360 WriteGPR_v100((* this), intermed); // this is for openbabel-1.100.2
361 istringstream interInput(intermed.str());
362
363 if (index == 0) // Automatic detection
364 {
365 if (!translator.CanExport(filename))
366 {
367 ErrorMessage(_("Cannot export that file type."));
368 return false;
369 }
370
371 ofile.open(filename, ios::out);
372 translator.Export(filename, interInput, ofile);
373 ofile.close();
374 }
375 else // By type picked by the user
376 {
377 ofile.open(filename, ios::out);
378 translator.Export(filename, index - 1, interInput, ofile);
379 ofile.close();
380 }
381
382 // change the original locale back...
383 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
384 setlocale(LC_NUMERIC, my_num_locale);
385
386 return true;
387 }
388
389 #endif // ENABLE_OPENBABEL
390
391 /*##############################################*/
392 /*##############################################*/
393
AddH(void)394 void project::AddH(void)
395 {
396 /* file_trans ft;
397 OBMol * obm = ft.CopyAll(this);
398 obm->AddHydrogens(false, false);
399 ft.Synchronize(); */
400
401 // above is the OpenBabel implementation of hydrogen adding.
402 // it seems to use bond length information to determine how many H's to add,
403 // which is problematic for protein X-ray structures (that often are not precise enough).
404 // TODO : make all these alternative add_h implementations available!!!!!!!!!!!!!!!
405
406 AddHydrogens(); // this is the library implementation...
407
408 ostringstream str;
409 str << _("Hydrogens added.") << endl << ends;
410
411 PrintToLog(str.str().c_str());
412 }
413
RemoveH(void)414 void project::RemoveH(void)
415 {
416 RemoveHydrogens();
417
418 ostringstream str;
419 str << _("Hydrogens removed.") << endl << ends;
420
421 PrintToLog(str.str().c_str());
422 }
423
GetDisplayListIDs(iGLu p1)424 iGLu project::GetDisplayListIDs(iGLu p1)
425 {
426 iGLu first = list_counter;
427 list_counter += p1;
428
429 return first;
430 }
431
DeleteDisplayLists(iGLu p1,iGLu p2)432 void project::DeleteDisplayLists(iGLu p1, iGLu p2)
433 {
434 for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
435 {
436 graphics_view_vector[n1]->GetWnd()->SetCurrent();
437 glDeleteLists(p1, p2);
438 }
439 }
440
441 /*##############################################*/
442 /*##############################################*/
443
AddAtom_lg(atom & p1)444 void project::AddAtom_lg(atom & p1)
445 {
446 custom_app * app2 = custom_app::GetAppC();
447
448 model::AddAtom_lg(p1);
449
450 atom_list.back().my_glname = app2->RegisterGLName(& atom_list.back());
451 app2->AtomAdded(& atom_list.back());
452 }
453
RemoveAtom(iter_al p1)454 void project::RemoveAtom(iter_al p1)
455 {
456 // first, discard ALL measure_tool information...
457
458 if (mt_a1 != NULL)
459 {
460 mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
461 mt_a1 = NULL;
462 }
463
464 if (mt_a2 != NULL)
465 {
466 mt_a2->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
467 mt_a2 = NULL;
468 }
469
470 if (mt_a3 != NULL)
471 {
472 mt_a3->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
473 mt_a3 = NULL;
474 }
475
476 // then proceed with the atom removal...
477
478 custom_app * app2 = custom_app::GetAppC();
479
480 app2->AtomRemoved(& (* p1));
481 app2->UnregisterGLNameByPtr(& (* p1));
482
483 model::RemoveAtom(p1);
484 }
485
AddBond(bond & p1)486 void project::AddBond(bond & p1)
487 {
488 custom_app * app2 = custom_app::GetAppC();
489
490 model::AddBond(p1);
491
492 app2->BondAdded(& bond_list.back());
493 }
494
RemoveBond(iter_bl p1)495 void project::RemoveBond(iter_bl p1)
496 {
497 custom_app * app2 = custom_app::GetAppC();
498
499 app2->BondRemoved(& (* p1));
500
501 model::RemoveBond(p1);
502 }
503
InvalidateGroups(void)504 void project::InvalidateGroups(void)
505 {
506 model::InvalidateGroups();
507 custom_app::GetAppC()->ClearChainsView();
508 }
509
UpdateChains(void)510 void project::UpdateChains(void)
511 {
512 model::UpdateChains();
513 custom_app::GetAppC()->BuildChainsView();
514 }
515
516 /*##############################################*/
517 /*##############################################*/
518
IsObject(const ogl_dummy_object * p1)519 i32s project::IsObject(const ogl_dummy_object * p1)
520 {
521 i32s index = NOT_DEFINED;
522 for (i32u n1 = 0;n1 < object_vector.size();n1++)
523 {
524 if (object_vector[n1] == p1) index = n1;
525 }
526
527 return index;
528 }
529
SelectObject(const ogl_dummy_object * p1)530 bool project::SelectObject(const ogl_dummy_object * p1)
531 {
532 i32s n1 = IsObject(p1);
533 if (n1 < 0) return false;
534
535 selected_object = object_vector[n1];
536
537 return true;
538 }
539
AddObject(ogl_smart_object * p1)540 void project::AddObject(ogl_smart_object * p1)
541 {
542 object_vector.push_back(p1);
543 selected_object = object_vector.back();
544
545 custom_app::GetAppC()->ObjectAdded(p1);
546 }
547
RemoveObject(ogl_dummy_object * p1)548 bool project::RemoveObject(ogl_dummy_object * p1)
549 {
550 i32s n1 = IsObject(p1);
551 if (n1 < 0) return false;
552
553 custom_app::GetAppC()->ObjectRemoved(object_vector[n1]);
554
555 object_vector.erase(object_vector.begin() + n1);
556 delete p1; return true;
557 }
558
559 // these are the measuring functions, that only take coordinates as input (so they are model-independent)...
560 // these are the measuring functions, that only take coordinates as input (so they are model-independent)...
561 // these are the measuring functions, that only take coordinates as input (so they are model-independent)...
562
measure_len(float * c1,float * c2)563 float measure_len(float * c1, float * c2)
564 {
565 v3d<float> v1(c1, c2);
566 return v1.len();
567 }
568
measure_ang(float * c1,float * c2,float * c3)569 float measure_ang(float * c1, float * c2, float * c3)
570 {
571 /* Vector v1, v2;
572 v1 = Vector(c1[0] - c2[0], c1[1] - c2[1], c1[2] - c2[2]);
573 v2 = Vector(c3[0] - c2[0], c3[1] - c2[1], c3[2] - c2[2]);
574 return VectorAngle(v1, v2); */
575
576 v3d<float> v1(c2, c1);
577 v3d<float> v2(c2, c3);
578 return v1.ang(v2) * 180.0 / M_PI;
579 }
580
measure_tor(float * c1,float * c2,float * c3,float * c4)581 float measure_tor(float * c1, float * c2, float * c3, float * c4)
582 {
583 /* Vector v1, v2, v3, v4;
584 v1 = Vector(c1[0], c1[1], c1[2]) * 10.0f;
585 v2 = Vector(c2[0], c2[1], c2[2]) * 10.0f;
586 v3 = Vector(c3[0], c3[1], c3[2]) * 10.0f;
587 v4 = Vector(c4[0], c4[1], c4[2]) * 10.0f;
588 return CalcTorsionAngle(v1, v2, v3, v4); */
589
590 v3d<float> v1(c2, c1);
591 v3d<float> v2(c2, c3);
592 v3d<float> v3(c3, c4);
593 return v1.tor(v2, v3) * 180.0 / M_PI;
594 }
595
596 /*##############################################*/
597 /*##############################################*/
598
AddGraphicsClient(custom_camera * cam,bool detached)599 oglview_wcl * project::AddGraphicsClient(custom_camera * cam, bool detached)
600 {
601 if (!cam)
602 {
603 fGL focus = GetDefaultFocus();
604 cam = new custom_camera(ogl_ol_static(), focus, this);
605
606 custom_app::GetAppC()->AddCamera(cam);
607
608 // also add a new light object by default...
609
610 ogl_light * l = new ogl_directional_light(ogl_ol_static());
611 custom_app::GetAppC()->AddLocalLight(l, cam);
612 }
613
614 oglview_wcl * wcl = new oglview_wcl(cam);
615
616 ostringstream title;
617 title << _("camera ") << cam->GetCCamI() << _(" window ") << wcl->my_wnd_number << ends;
618 wcl->SetTitle(title.str().c_str());
619
620 base_wnd * wnd = CreateGraphicsWnd(detached);
621
622 wcl->LinkWnd(wnd);
623
624 graphics_view_vector.push_back(wcl);
625 custom_app::GetAppC()->GraphicsClientAdded(wcl);
626
627 custom_app::GetAppC()->SetupLights(cam);
628 custom_app::GetAppC()->UpdateAllWindowTitles();
629
630 return wcl;
631 }
632
RemoveGraphicsClient(oglview_wcl * wcl,bool force)633 bool project::RemoveGraphicsClient(oglview_wcl * wcl, bool force)
634 {
635 ogl_camera * cam = wcl->GetCam();
636
637 i32s views_for_other_cams = 0;
638 i32s other_views_for_this_cam = 0;
639
640 i32s index = NOT_DEFINED;
641 for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
642 {
643 if (graphics_view_vector[n1] == wcl)
644 {
645 index = n1;
646 continue;
647 }
648
649 if (graphics_view_vector[n1]->GetCam() != cam)
650 {
651 views_for_other_cams++;
652 }
653 else
654 {
655 other_views_for_this_cam++;
656 }
657 }
658
659 if (index < 0) assertion_failed(__FILE__, __LINE__, "index < 0");
660
661 if (!force && (views_for_other_cams + other_views_for_this_cam < 1)) // refuse to close the last view!!!
662 {
663 ErrorMessage(_("This is the last graphics view for\nthis project - can't close it."));
664 return false;
665 }
666
667 // now let's remove the window and the client...
668
669 base_wnd * wnd = wcl->GetWnd();
670
671 wcl->UnlinkWnd();
672
673 DestroyGraphicsWnd(wnd);
674 wnd = NULL;
675
676 custom_app::GetAppC()->GraphicsClientRemoved(wcl);
677
678 graphics_view_vector.erase(graphics_view_vector.begin() + index);
679
680 delete wcl;
681 wcl = NULL;
682
683 // now we can also remove the camera, if needed...
684
685 if (!other_views_for_this_cam)
686 {
687 custom_app::GetAppC()->RemoveCamera(cam);
688
689 delete cam;
690 cam = NULL;
691
692 // disable selected_object since it could have been
693 // a light object that just got deleted ; FIX_ME_LATER
694
695 selected_object = NULL; // ???
696 }
697
698 custom_app::GetAppC()->UpdateAllWindowTitles();
699
700 return true;
701 }
702
IsThisLastGraphicsClient(oglview_wcl * wcl)703 bool project::IsThisLastGraphicsClient(oglview_wcl * wcl)
704 {
705 ogl_camera * cam = wcl->GetCam();
706
707 i32s views_for_other_cams = 0;
708 i32s other_views_for_this_cam = 0;
709
710 for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
711 {
712 if (graphics_view_vector[n1] == wcl) continue;
713
714 if (graphics_view_vector[n1]->GetCam() != cam)
715 {
716 views_for_other_cams++;
717 }
718 else
719 {
720 other_views_for_this_cam++;
721 }
722 }
723
724 if (views_for_other_cams + other_views_for_this_cam < 1) return true;
725 else return false;
726 }
727
AddPlot1DClient(const char * s1,const char * sv,bool detached)728 p1dview_wcl * project::AddPlot1DClient(const char * s1, const char * sv, bool detached)
729 {
730 p1dview_wcl * wcl = new p1dview_wcl(s1, sv);
731 base_wnd * wnd = CreatePlot1DWnd(detached);
732
733 wcl->LinkWnd(wnd);
734
735 plotting_view_vector.push_back(wcl);
736 custom_app::GetAppC()->PlottingClientAdded(wcl);
737
738 custom_app::GetAppC()->UpdateAllWindowTitles();
739
740 return wcl;
741 }
742
AddPlot2DClient(const char * s1,const char * s2,const char * sv,bool detached)743 p2dview_wcl * project::AddPlot2DClient(const char * s1, const char * s2, const char * sv, bool detached)
744 {
745 p2dview_wcl * wcl = new p2dview_wcl(s1, s2, sv);
746 base_wnd * wnd = CreatePlot2DWnd(detached);
747
748 wcl->LinkWnd(wnd);
749
750 plotting_view_vector.push_back(wcl);
751 custom_app::GetAppC()->PlottingClientAdded(wcl);
752
753 custom_app::GetAppC()->UpdateAllWindowTitles();
754
755 return wcl;
756 }
757
AddEnergyLevelDiagramClient(bool detached)758 eldview_wcl * project::AddEnergyLevelDiagramClient(bool detached)
759 {
760 eldview_wcl * wcl = new eldview_wcl();
761 base_wnd * wnd = CreateEnergyLevelDiagramWnd(detached);
762
763 wcl->LinkWnd(wnd);
764
765 plotting_view_vector.push_back(wcl);
766 custom_app::GetAppC()->PlottingClientAdded(wcl);
767
768 custom_app::GetAppC()->UpdateAllWindowTitles();
769
770 return wcl;
771 }
772
AddReactionCoordinatePlotClient(const char * s1,const char * sv,bool detached)773 rcpview_wcl * project::AddReactionCoordinatePlotClient(const char * s1, const char * sv, bool detached)
774 {
775 rcpview_wcl * wcl = new rcpview_wcl(s1, sv);
776 base_wnd * wnd = CreateReactionCoordinatePlotWnd(detached);
777
778 wcl->LinkWnd(wnd);
779
780 plotting_view_vector.push_back(wcl);
781 custom_app::GetAppC()->PlottingClientAdded(wcl);
782
783 custom_app::GetAppC()->UpdateAllWindowTitles();
784
785 return wcl;
786 }
787
AddGenericProteinChainClient(bool detached)788 gpcview_wcl * project::AddGenericProteinChainClient(bool detached)
789 {
790 gpcview_wcl * wcl = new gpcview_wcl();
791 base_wnd * wnd = CreateGenericProteinChainWnd(detached);
792
793 wcl->LinkWnd(wnd);
794
795 plotting_view_vector.push_back(wcl);
796 custom_app::GetAppC()->PlottingClientAdded(wcl);
797
798 custom_app::GetAppC()->UpdateAllWindowTitles();
799
800 return wcl;
801 }
802
RemovePlottingClient(base_wcl * wcl)803 bool project::RemovePlottingClient(base_wcl * wcl)
804 {
805 i32s index = NOT_DEFINED;
806 for (i32u n2 = 0;n2 < plotting_view_vector.size();n2++)
807 {
808 if (plotting_view_vector[n2] == wcl) index = n2;
809 }
810
811 if (index < 0) assertion_failed(__FILE__, __LINE__, "index < 0");
812
813 // now let's remove the window and the client...
814
815 base_wnd * wnd = wcl->GetWnd();
816
817 wcl->UnlinkWnd();
818
819 DestroyPlottingWnd(wnd);
820 wnd = NULL;
821
822 custom_app::GetAppC()->PlottingClientRemoved(wcl);
823
824 plotting_view_vector.erase(plotting_view_vector.begin() + index);
825
826 delete wcl;
827 wcl = NULL;
828
829 return true;
830 }
831
UpdateAllViews(void)832 void project::UpdateAllViews(void)
833 {
834 UpdateAllGraphicsViews();
835 UpdateAllPlottingViews();
836
837 // the project view, if exists, will take
838 // care of itself and does not need updating.
839 }
840
UpdateAllGraphicsViews(bool flag)841 void project::UpdateAllGraphicsViews(bool flag)
842 {
843 for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
844 {
845 graphics_view_vector[n1]->GetWnd()->RequestUpdate(flag);
846 }
847 }
848
UpdateAllPlottingViews(bool flag)849 void project::UpdateAllPlottingViews(bool flag)
850 {
851 for (i32u n1 = 0;n1 < plotting_view_vector.size();n1++)
852 {
853 plotting_view_vector[n1]->GetWnd()->RequestUpdate(flag);
854 }
855 }
856
UpdateGraphicsViews(ogl_camera * cam,bool flag)857 void project::UpdateGraphicsViews(ogl_camera * cam, bool flag)
858 {
859 for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
860 {
861 if (graphics_view_vector[n1]->GetCam() != cam) continue;
862 graphics_view_vector[n1]->GetWnd()->RequestUpdate(flag);
863 }
864 }
865
UpdateGraphicsView(oglview_wcl * wcl,bool flag)866 void project::UpdateGraphicsView(oglview_wcl * wcl, bool flag)
867 {
868 wcl->GetWnd()->RequestUpdate(flag);
869 }
870
871 /*##############################################*/
872 /*##############################################*/
873
ProcessCommandString(oglview_wcl * wcl,const char * command)874 void project::ProcessCommandString(oglview_wcl * wcl, const char * command)
875 {
876 ostringstream str1;
877 str1 << _("Processing Command : ") << command << endl << ends;
878 PrintToLog(str1.str().c_str());
879
880 istringstream istr(command);
881 char kw1[32]; istr >> kw1; // the 1st keyword.
882
883 if (!strcmp("help", kw1))
884 {
885 ostringstream str;
886
887 str << _("> AVAILABLE COMMANDS:") << endl; // use alphabetical order???
888
889 str << _("> add light (local/global) (directional/spotlight) -- add a new light object.") << endl;
890
891 str << _("> add plane <vf> <cf> <cscale1> <AUTO/cscale2> <dim> <res> <tp> <alpha> -- add a plane object.") << endl;
892 str << _("> where: <vf> = value function : esp vdws eldens mo mod unity") << endl;
893 str << _("> <cf> = colour function : red green blue rb1 rb2") << endl;
894 str << _("> <cscale1> = scaling value for calculating the colours") << endl;
895 str << _("> <cscale2> = scaling offset for calculating the colours") << endl;
896 str << _("> <dim> = dimension of the plane object (in nm units)") << endl;
897 str << _("> <res> = resolution of the plane object") << endl;
898 str << _("> <tp> = 0 or 1 telling if the object is transparent") << endl;
899 str << _("> <alpha> = transparency alpha value") << endl;
900
901 str << _("> add volrend <vf> <cf> <cscale1> <AUTO/cscale2> <dim> <res> <alpha> -- add a volume-rendering object.") << endl;
902 str << _("> where: <vf> = value function : esp vdws eldens mo mod unity") << endl;
903 str << _("> <cf> = colour function : red green blue rb1 rb2") << endl;
904 str << _("> <cscale1> = scaling value for calculating the colours") << endl;
905 str << _("> <cscale2> = scaling offset for calculating the colours") << endl;
906 str << _("> <dim> = dimension of the plane object (in nm units)") << endl;
907 str << _("> <res> = resolution of the plane object") << endl;
908 str << _("> <alpha> = transparency alpha value") << endl;
909
910 str << _("> add surf1 <vf1> <vf2> <cf> <sscale> <cscale1> <AUTO/cscale2> <dim> <res> <solid> <tp> <alpha> -- add a single surface object.") << endl;
911 str << _("> where: <vf1> = value function for calculating the surface : esp vdws eldens mo mod unity") << endl;
912 str << _("> <vf2> = value function for calculating the colours : esp vdws eldens mo mod unity") << endl;
913 str << _("> <cf> = colour function : red green blue rb1 rb2") << endl;
914 str << _("> <sscale> = scaling value for calculating the surface") << endl;
915 str << _("> <cscale1> = scaling value for calculating the colours") << endl;
916 str << _("> <cscale2> = scaling offset for calculating the colours") << endl;
917 str << _("> <dim> = dimension of the plane object (in nm units)")<< endl;
918 str << _("> <res> = resolution of the plane object") << endl;
919 str << _("> <solid> = 0 or 1 telling if the object is solid") << endl;
920 str << _("> <tp> = 0 or 1 telling if the object is transparent")<< endl;
921 str << _("> <alpha> = transparency alpha value") << endl;
922
923 str << _("> add surf2 <vf1> <vf2> <cf1> <cf2> <sscale1> <sscale2> <cscale1> <AUTO/cscale2> <dim> <res> <solid> <tp> <alpha> -- add a pair of surface objects.") << endl;
924 str << _("> where: <vf1> = value function for calculating the surface : esp vdws eldens mo mod unity") << endl;
925 str << _("> <vf2> = value function for calculating the colours : esp vdws eldens mo mod unity") << endl;
926 str << _("> <cf1> = colour function for 1st surface : red green blue rb1 rb2") << endl;
927 str << _("> <cf2> = colour function for 2nd surface : red green blue rb1 rb2") << endl;
928 str << _("> <sscale1> = scaling value for calculating the surface for 1st surface") << endl;
929 str << _("> <sscale2> = scaling value for calculating the surface for 2nd surface")<< endl;
930 str << _("> <cscale1> = scaling value for calculating the colours")<< endl;
931 str << _("> <cscale2> = scaling offset for calculating the colours")<< endl;
932 str << _("> <dim> = dimension of the plane object (in nm units)")<< endl;
933 str << _("> <res> = resolution of the plane object")<< endl;
934 str << _("> <solid> = 0 or 1 telling if the object is solid")<< endl;
935 str << _("> <tp> = 0 or 1 telling if the object is transparent")<< endl;
936 str << _("> <alpha> = transparency alpha value") << endl;
937
938 str << _("> help -- print all available commands in command strings.") << endl;
939
940 str << _("> energy -- calculate a single-point energy.") << endl;
941 str << _("> geom_opt -- do a geometry optimization run using default options.") << endl;
942 str << _("> mol_dyn -- do a molecular dynamics run using default options.") << endl;
943
944 str << _("> random_search <cycles> <optsteps> -- perform a random conformational search.") << endl;
945 str << _("> systematic_search <divisions> <optsteps> -- perform a systematic conformational search.") << endl;
946 str << _("> montecarlo_search <init_cycles> <simul_cycles> <optsteps> -- perform a MonteCarlo search.") << endl;
947
948 str << _("> make_plot1 A B C D <div> <start_ang> <end_ang> <optsteps> -- create a 1D energy vs. torsion plot.") << endl;
949 str << _("> make_plot2 A B C D <div> <start_ang> <end_ang> I J K L <div> <start_ang> <end_ang> <optsteps> -- create a 2D energy vs. torsions plot.") << endl;
950
951 str << _("> population_analysis_ESP -- determine atomic charges using an ESP fit (for QM methods only).") << endl;
952
953 str << _("> transition_state_search <delta_e> <initial_fc> -- perform a transition state search (for QM methods only).") << endl;
954 str << _("> stationary_state_search <steps> -- perform a search for a structure with no forces.") << endl;
955
956 str << _("> set_current_orbital <orbital_index> -- set the current orbtal index for plotting the orbitals.") << endl;
957
958 str << _("> update_chains -- detect polymer chains using sequence builder.") << endl;
959 str << _("> build_amino <sequence> (helix/strand) -- amino acid sequence builder.") << endl;
960 str << _("> build_nucleic <sequence> -- nucleic acid sequence builder.") << endl;
961
962 str << _("> orient <crdset> -- orient the system in the XYZ coordinate system.") << endl;
963
964 str << _("> solvate_box <x-hdim> <y-hdim> <z-hdim> (<density> <filename> (export)) -- setup a solvation box.") << endl;
965 str << _("> solvate_sphere <rad_solute> <rad_solvent> (<density> <filename>) -- setup a solvation sphere.")<< endl;
966
967 str << _("> set_formal_charge <index> <charge> -- set formal charges to atoms.") << endl;
968
969 str << _("> evaluate_Bfact -- evaluate B-factors for selected atoms (a trajectory file must be open).") << endl;
970 str << _("> evaluate_diffconst <dt> -- evaluate diffusion constants for selected atoms (a trajectory file must be open, dt = time difference between frames [fs]).") << endl;
971
972 str << ends;
973 PrintToLog(str.str().c_str());
974 return;
975 }
976
977 if (!strcmp("add", kw1))
978 {
979 char kw2[32]; istr >> kw2; // the 2nd keyword; type of the object to add.
980
981 if (!strcmp("light", kw2))
982 {
983 char kw3[32]; istr >> kw3;
984 char kw4[32]; istr >> kw4;
985
986 bool is_local = true;
987 bool is_directional = true;
988
989 if (kw3[0] == 'g' || kw4[0] == 'g') is_local = false; // global
990 if (kw3[0] == 's' || kw4[0] == 's') is_directional = false; // spotlight
991
992 ogl_light * new_light;
993 if (is_directional) new_light = new ogl_directional_light(ogl_ol_static());
994 else
995 {
996 new_light = new rendered_spot_light(ogl_ol_static());
997
998 const fGL trans[3] = { 0.0, 0.0, -1.0 };
999 new_light->TranslateObject((const fGL *) trans, wcl->GetCam()->GetSafeLD());
1000 }
1001
1002 base_app * app = base_app::GetAppB();
1003 if (!is_local) app->AddGlobalLight(new_light);
1004 else app->AddLocalLight(new_light, wcl->GetCam());
1005
1006 if (!is_local || oglview_wcl::draw_info) UpdateAllGraphicsViews();
1007 else if (is_local) UpdateGraphicsViews(wcl->GetCam());
1008
1009 ostringstream strR;
1010 strR << _("Added a new object : light (");
1011 strR << (is_local ? _("local") : _("global")) << " ";
1012 strR << (is_directional ? _("directional") : _("spotlight")) << ")." << endl << ends;
1013 PrintToLog(strR.str().c_str());
1014 return;
1015 }
1016
1017 if (!strcmp("plane", kw2))
1018 {
1019 char kw3[32]; istr >> kw3;
1020 char kw4[32]; istr >> kw4;
1021 char kw5[32]; istr >> kw5;
1022 char kw6[32]; istr >> kw6;
1023 char kw7[32]; istr >> kw7;
1024 char kw8[32]; istr >> kw8;
1025 char kw9[32]; istr >> kw9;
1026 char kwA[32]; istr >> kwA;
1027 char ** endptr = NULL;
1028
1029 ogl_cp_param cpp;
1030 cpp.prj = this; cpp.ref = GetCurrentSetup()->GetCurrentEngine();
1031
1032 if (!strcmp(kw3, "esp")) cpp.vf = (ValueFunction *) value_ESP;
1033 else if (!strcmp(kw3, "vdws")) cpp.vf = (ValueFunction *) value_VDWSurf;
1034 else if (!strcmp(kw3, "eldens")) cpp.vf = (ValueFunction *) value_ElDens;
1035 else if (!strcmp(kw3, "mo")) cpp.vf = (ValueFunction *) value_Orbital;
1036 else if (!strcmp(kw3, "mod")) cpp.vf = (ValueFunction *) value_OrbDens;
1037 else if (!strcmp(kw3, "unity")) cpp.vf = (ValueFunction *) GetUnity;
1038 else
1039 {
1040 ostringstream strE;
1041 strE << _("ERROR : add plane : unknown value function ") << kw3 << "." << endl << ends;
1042
1043 PrintToLog(strE.str().c_str());
1044 return;
1045 }
1046
1047 if (!strcmp(kw4, "red")) cpp.cf = (ColorFunction *) GetRedColor;
1048 else if (!strcmp(kw4, "green")) cpp.cf = (ColorFunction *) GetGreenColor;
1049 else if (!strcmp(kw4, "blue")) cpp.cf = (ColorFunction *) GetBlueColor;
1050 else if (!strcmp(kw4, "rb1")) cpp.cf = (ColorFunction *) GetRBRange1;
1051 else if (!strcmp(kw4, "rb2")) cpp.cf = (ColorFunction *) GetRBRange2;
1052 else
1053 {
1054 ostringstream strE;
1055 strE << _("ERROR : add plane : unknown colour function ") << kw4 << "." << endl << ends;
1056
1057 PrintToLog(strE.str().c_str());
1058 return;
1059 }
1060
1061 f64 cscale1 = strtod(kw5, endptr);
1062
1063 f64 cscale2 = 0.0; bool auto_cv2 = false;
1064 if (!strcmp(kw6, "AUTO")) auto_cv2 = true;
1065 else cscale2 = strtod(kw6, endptr);
1066
1067 f64 dim = strtod(kw7, endptr);
1068
1069 i32s res = strtol(kw8, endptr, 10);
1070 if (res < 2) res = 2;
1071
1072 i32s tp = strtol(kw9, endptr, 10);
1073 if (tp < 0) tp = 0; if (tp > 1) tp = 1;
1074
1075 f64 alpha = strtod(kwA, endptr);
1076
1077 cpp.dim = dim; cpp.np = res;
1078 cpp.transparent = tp; cpp.automatic_cv2 = auto_cv2;
1079
1080 cpp.cvalue1 = cscale1;
1081 cpp.cvalue2 = cscale2;
1082 cpp.alpha = alpha;
1083
1084 ostringstream strN;
1085 strN << kw3 << "-" << ends;
1086
1087 AddObject(new ogl_color_plane_object(ogl_ol_static(), cpp, strN.str().c_str()));
1088 UpdateAllGraphicsViews();
1089
1090 ostringstream strR;
1091 strR << _("Added a new object : plane (") << kw3 << " " << kw4 << ")." << endl << ends;
1092 PrintToLog(strR.str().c_str());
1093 return;
1094 }
1095
1096 if (!strcmp("volrend", kw2))
1097 {
1098 char kw3[32]; istr >> kw3;
1099 char kw4[32]; istr >> kw4;
1100 char kw5[32]; istr >> kw5;
1101 char kw6[32]; istr >> kw6;
1102 char kw7[32]; istr >> kw7;
1103 char kw8[32]; istr >> kw8;
1104 char kw9[32]; istr >> kw9;
1105 char ** endptr = NULL;
1106
1107 ogl_cp_param cpp;
1108 cpp.prj = this; cpp.ref = GetCurrentSetup()->GetCurrentEngine();
1109
1110 if (!strcmp(kw3, "esp")) cpp.vf = (ValueFunction *) value_ESP;
1111 else if (!strcmp(kw3, "vdws")) cpp.vf = (ValueFunction *) value_VDWSurf;
1112 else if (!strcmp(kw3, "eldens")) cpp.vf = (ValueFunction *) value_ElDens;
1113 else if (!strcmp(kw3, "mo")) cpp.vf = (ValueFunction *) value_Orbital;
1114 else if (!strcmp(kw3, "mod")) cpp.vf = (ValueFunction *) value_OrbDens;
1115 else if (!strcmp(kw3, "unity")) cpp.vf = (ValueFunction *) GetUnity;
1116 else
1117 {
1118 ostringstream strE;
1119 strE << _("ERROR : add volrend : unknown value function ") << kw3 << "." << endl << ends;
1120
1121 PrintToLog(strE.str().c_str());
1122 return;
1123 }
1124
1125 if (!strcmp(kw4, "red")) cpp.cf = (ColorFunction *) GetRedColor;
1126 else if (!strcmp(kw4, "green")) cpp.cf = (ColorFunction *) GetGreenColor;
1127 else if (!strcmp(kw4, "blue")) cpp.cf = (ColorFunction *) GetBlueColor;
1128 else if (!strcmp(kw4, "rb1")) cpp.cf = (ColorFunction *) GetRBRange1;
1129 else if (!strcmp(kw4, "rb2")) cpp.cf = (ColorFunction *) GetRBRange2;
1130 else
1131 {
1132 ostringstream strE;
1133 strE << _("ERROR : add volrend : unknown colour function ") << kw4 <<"." << endl << ends;
1134
1135 PrintToLog(strE.str().c_str());
1136 return;
1137 }
1138
1139 f64 cscale1 = strtod(kw5, endptr);
1140
1141 f64 cscale2 = 0.0; bool auto_cv2 = false;
1142 if (!strcmp(kw6, "AUTO")) auto_cv2 = true;
1143 else cscale2 = strtod(kw6, endptr);
1144
1145 f64 dim = strtod(kw7, endptr);
1146
1147 i32s res = strtol(kw8, endptr, 10);
1148 if (res < 4) res = 4;
1149
1150 f64 alpha = strtod(kw9, endptr);
1151
1152 cpp.dim = dim; cpp.np = res;
1153 cpp.transparent = true; cpp.automatic_cv2 = auto_cv2;
1154
1155 cpp.cvalue1 = cscale1;
1156 cpp.cvalue2 = cscale2;
1157 cpp.alpha = alpha;
1158
1159 ostringstream strN;
1160 strN << kw3 << "-" << ends;
1161
1162 AddObject(new ogl_volume_rendering_object(ogl_ol_static(), cpp, res / 2, dim / 2.0, (* wcl->GetCam()), strN.str().c_str()));
1163 UpdateAllGraphicsViews();
1164
1165 ostringstream strR;
1166 strR << _("Added a new object : volrend (") << kw3 << " " << kw4 << ")." << endl << ends;
1167 PrintToLog(strR.str().c_str());
1168 return;
1169 }
1170
1171 if (!strcmp("surf1", kw2))
1172 {
1173 char kw3[32]; istr >> kw3; // vf1
1174 char kw4[32]; istr >> kw4; // vf2
1175 char kw5[32]; istr >> kw5; // cf
1176 char kw6[32]; istr >> kw6; // sscale
1177 char kw7[32]; istr >> kw7; // cscale1
1178 char kw8[32]; istr >> kw8; // AUTO/cscale2
1179 char kw9[32]; istr >> kw9; // dim
1180 char kwA[32]; istr >> kwA; // res
1181 char kwB[32]; istr >> kwB; // solid
1182 char kwC[32]; istr >> kwC; // tp
1183 char kwD[32]; istr >> kwD; // alpha
1184 char ** endptr = NULL;
1185
1186 ogl_cs_param csp1;
1187
1188 csp1.prj = this; csp1.ref = GetCurrentSetup()->GetCurrentEngine(); csp1.next = NULL;
1189
1190 if (!strcmp(kw3, "esp")) csp1.vf1 = (ValueFunction *) value_ESP;
1191 else if (!strcmp(kw3, "vdws")) csp1.vf1 = (ValueFunction *) value_VDWSurf;
1192 else if (!strcmp(kw3, "eldens")) csp1.vf1 = (ValueFunction *) value_ElDens;
1193 else if (!strcmp(kw3, "mo")) csp1.vf1 = (ValueFunction *) value_Orbital;
1194 else if (!strcmp(kw3, "mod")) csp1.vf1 = (ValueFunction *) value_OrbDens;
1195 else if (!strcmp(kw3, "unity")) csp1.vf1 = (ValueFunction *) GetUnity;
1196 else
1197 {
1198 ostringstream strE;
1199 strE << _("ERROR : add surf1 : unknown value function 1 ") << kw3 << "." << endl << ends;
1200
1201 PrintToLog(strE.str().c_str());
1202 return;
1203 }
1204
1205 if (!strcmp(kw4, "esp")) csp1.vf2 = (ValueFunction *) value_ESP;
1206 else if (!strcmp(kw4, "vdws")) csp1.vf2 = (ValueFunction *) value_VDWSurf;
1207 else if (!strcmp(kw4, "eldens")) csp1.vf2 = (ValueFunction *) value_ElDens;
1208 else if (!strcmp(kw4, "mo")) csp1.vf2 = (ValueFunction *) value_Orbital;
1209 else if (!strcmp(kw4, "mod")) csp1.vf2 = (ValueFunction *) value_OrbDens;
1210 else if (!strcmp(kw4, "unity")) csp1.vf2 = (ValueFunction *) GetUnity;
1211 else
1212 {
1213 ostringstream strE;
1214 strE << _("ERROR : add surf1 : unknown value function 2 ") << kw4 << "." << endl << ends;
1215
1216 PrintToLog(strE.str().c_str());
1217 return;
1218 }
1219
1220 if (!strcmp(kw5, "red")) csp1.cf = (ColorFunction *) GetRedColor;
1221 else if (!strcmp(kw5, "green")) csp1.cf = (ColorFunction *) GetGreenColor;
1222 else if (!strcmp(kw5, "blue")) csp1.cf = (ColorFunction *) GetBlueColor;
1223 else if (!strcmp(kw5, "rb1")) csp1.cf = (ColorFunction *) GetRBRange1;
1224 else if (!strcmp(kw5, "rb2")) csp1.cf = (ColorFunction *) GetRBRange2;
1225 else
1226 {
1227 ostringstream strE;
1228 strE << _("ERROR : add surf1 : unknown colour function ") << kw5 << "." << endl << ends;
1229
1230 PrintToLog(strE.str().c_str());
1231 return;
1232 }
1233
1234 f64 sscale = strtod(kw6, endptr);
1235 f64 cscale1 = strtod(kw7, endptr);
1236
1237 f64 cscale2 = 0.0; bool auto_cv2 = false;
1238 if (!strcmp(kw8, "AUTO")) auto_cv2 = true;
1239 else cscale2 = strtod(kw8, endptr);
1240
1241 f64 dim = strtod(kw9, endptr);
1242
1243 i32s res = strtol(kwA, endptr, 10);
1244 if (res < 4) res = 4;
1245
1246 i32s solid = strtol(kwB, endptr, 10);
1247 if (solid < 0) solid = 0; if (solid > 1) solid = 1;
1248
1249 i32s tp = strtol(kwC, endptr, 10);
1250 if (tp < 0) tp = 0; if (tp > 1) tp = 1;
1251
1252 f64 alpha = strtod(kwD, endptr);
1253
1254 static fGL dim_arr[3];
1255 dim_arr[0] = dim_arr[1] = dim_arr[2] = dim;
1256
1257 static i32s res_arr[3];
1258 res_arr[0] = res_arr[1] = res_arr[2] = res;
1259
1260 csp1.dim = dim_arr; csp1.np = res_arr;
1261 csp1.transparent = tp; csp1.automatic_cv2 = auto_cv2; csp1.wireframe = !solid;
1262
1263 csp1.svalue = sscale;
1264 csp1.cvalue1 = cscale1;
1265 csp1.cvalue2 = cscale2;
1266 csp1.alpha = alpha;
1267
1268 csp1.toler = fabs(1.0e-6 * sscale); csp1.maxc = 250;
1269
1270 ostringstream strN;
1271 strN << kw3 << "-" << kw4 << "-" << ends;
1272
1273 AddObject(new ogl_color_surface_object(ogl_ol_static(), csp1, strN.str().c_str()));
1274 UpdateAllGraphicsViews();
1275
1276 ostringstream strR;
1277 strR << _("Added a new object : surf1 (") << kw3 << " " << kw4 << " " << kw5 << ")." << endl << ends;
1278 PrintToLog(strR.str().c_str());
1279 return;
1280 }
1281
1282 if (!strcmp("surf2", kw2))
1283 {
1284 char kw3[32]; istr >> kw3; // vf1
1285 char kw4[32]; istr >> kw4; // vf2
1286 char kw5[32]; istr >> kw5; // cf1
1287 char kw6[32]; istr >> kw6; // cf2
1288 char kw7[32]; istr >> kw7; // sscale1
1289 char kw8[32]; istr >> kw8; // sscale2
1290 char kw9[32]; istr >> kw9; // cscale1
1291 char kwA[32]; istr >> kwA; // AUTO/cscale2
1292 char kwB[32]; istr >> kwB; // dim
1293 char kwC[32]; istr >> kwC; // res
1294 char kwD[32]; istr >> kwD; // solid
1295 char kwE[32]; istr >> kwE; // tp
1296 char kwF[32]; istr >> kwF; // alpha
1297 char ** endptr = NULL;
1298
1299 ogl_cs_param csp2a; ogl_cs_param csp2b;
1300
1301 csp2a.prj = this; csp2a.ref = GetCurrentSetup()->GetCurrentEngine(); csp2a.next = & csp2b;
1302
1303 csp2b.prj = this; csp2b.ref = GetCurrentSetup()->GetCurrentEngine(); csp2b.next = NULL;
1304
1305 if (!strcmp(kw3, "esp")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_ESP;
1306 else if (!strcmp(kw3, "vdws")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_VDWSurf;
1307 else if (!strcmp(kw3, "eldens")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_ElDens;
1308 else if (!strcmp(kw3, "mo")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_Orbital;
1309 else if (!strcmp(kw3, "mod")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_OrbDens;
1310 else if (!strcmp(kw3, "unity")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) GetUnity;
1311 else
1312 {
1313 ostringstream strE;
1314 strE << _("ERROR : add surf2 : unknown value function 1 ") << kw3 << "." << endl << ends;
1315
1316 PrintToLog(strE.str().c_str());
1317 return;
1318 }
1319
1320 if (!strcmp(kw4, "esp")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_ESP;
1321 else if (!strcmp(kw4, "vdws")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_VDWSurf;
1322 else if (!strcmp(kw4, "eldens")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_ElDens;
1323 else if (!strcmp(kw4, "mo")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_Orbital;
1324 else if (!strcmp(kw4, "mod")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_OrbDens;
1325 else if (!strcmp(kw4, "unity")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) GetUnity;
1326 else
1327 {
1328 ostringstream strE;
1329 strE << _("ERROR : add surf2 : unknown value function 2 ") << kw4 << "." << endl << ends;
1330
1331 PrintToLog(strE.str().c_str());
1332 return;
1333 }
1334
1335 if (!strcmp(kw5, "red")) csp2a.cf = (ColorFunction *) GetRedColor;
1336 else if (!strcmp(kw5, "green")) csp2a.cf = (ColorFunction *) GetGreenColor;
1337 else if (!strcmp(kw5, "blue")) csp2a.cf = (ColorFunction *) GetBlueColor;
1338 else if (!strcmp(kw5, "rb1")) csp2a.cf = (ColorFunction *) GetRBRange1;
1339 else if (!strcmp(kw5, "rb2")) csp2a.cf = (ColorFunction *) GetRBRange2;
1340 else
1341 {
1342 ostringstream strE;
1343 strE << _("ERROR : add surf2 : unknown colour function 1 ") << kw5 << "." << endl << ends;
1344
1345 PrintToLog(strE.str().c_str());
1346 return;
1347 }
1348
1349 if (!strcmp(kw6, "red")) csp2b.cf = (ColorFunction *) GetRedColor;
1350 else if (!strcmp(kw6, "green")) csp2b.cf = (ColorFunction *) GetGreenColor;
1351 else if (!strcmp(kw6, "blue")) csp2b.cf = (ColorFunction *) GetBlueColor;
1352 else if (!strcmp(kw6, "rb1")) csp2b.cf = (ColorFunction *) GetRBRange1;
1353 else if (!strcmp(kw6, "rb2")) csp2b.cf = (ColorFunction *) GetRBRange2;
1354 else
1355 {
1356 ostringstream strE;
1357 strE << _("ERROR : add surf2 : unknown colour function 2 ") << kw6 << "." << endl << ends;
1358
1359 PrintToLog(strE.str().c_str());
1360 return;
1361 }
1362
1363 f64 sscale1 = strtod(kw7, endptr);
1364 f64 sscale2 = strtod(kw8, endptr);
1365 f64 cscale1 = strtod(kw9, endptr);
1366
1367 f64 cscale2 = 0.0; bool auto_cv2 = false;
1368 if (!strcmp(kwA, "AUTO")) auto_cv2 = true;
1369 else cscale2 = strtod(kwA, endptr);
1370
1371 f64 dim = strtod(kwB, endptr);
1372
1373 i32s res = strtol(kwC, endptr, 10);
1374 if (res < 4) res = 4;
1375
1376 i32s solid = strtol(kwD, endptr, 10);
1377 if (solid < 0) solid = 0; if (solid > 1) solid = 1;
1378
1379 i32s tp = strtol(kwE, endptr, 10);
1380 if (tp < 0) tp = 0; if (tp > 1) tp = 1;
1381
1382 f64 alpha = strtod(kwF, endptr);
1383
1384 static fGL dim_arr[3];
1385 dim_arr[0] = dim_arr[1] = dim_arr[2] = dim;
1386
1387 static i32s res_arr[3];
1388 res_arr[0] = res_arr[1] = res_arr[2] = res;
1389
1390 csp2a.dim = dim_arr; csp2a.np = res_arr;
1391 csp2a.transparent = tp; csp2a.automatic_cv2 = auto_cv2; csp2a.wireframe = !solid;
1392
1393 csp2a.svalue = sscale1;
1394 csp2a.cvalue1 = cscale1;
1395 csp2a.cvalue2 = cscale2;
1396 csp2a.alpha = alpha;
1397
1398 csp2a.toler = fabs(1.0e-6 * sscale1); csp2a.maxc = 250;
1399
1400 csp2b.dim = dim_arr; csp2b.np = res_arr;
1401 csp2b.transparent = tp; csp2b.automatic_cv2 = auto_cv2; csp2b.wireframe = !solid;
1402
1403 csp2b.svalue = sscale2;
1404 csp2b.cvalue1 = cscale1;
1405 csp2b.cvalue2 = cscale2;
1406 csp2b.alpha = alpha;
1407
1408 csp2b.toler = fabs(1.0e-6 * sscale2); csp2b.maxc = 250;
1409
1410 ostringstream strN;
1411 strN << kw3 << "-" << kw4 << "-" << ends;
1412
1413 AddObject(new ogl_color_surface_object(ogl_ol_static(), csp2a, strN.str().c_str()));
1414 UpdateAllGraphicsViews();
1415
1416 ostringstream strR;
1417 strR << _("Added a new object : surf2 (") << kw3 << " " << kw4 << " " << kw5 << " " << kw6 << ")." << endl << ends;
1418 PrintToLog(strR.str().c_str());
1419 return;
1420 }
1421
1422 ostringstream strE;
1423 strE << _("ERROR : could not process command \"add\" for parameter ") << kw2 << "." << endl << ends;
1424
1425 PrintToLog(strE.str().c_str());
1426 return;
1427 }
1428
1429 if (!strcmp("energy", kw1))
1430 {
1431 DoEnergy();
1432 return;
1433 }
1434
1435 if (!strcmp("geom_opt", kw1)) // todo: how to set the options here?
1436 {
1437 setup * su = GetCurrentSetup();
1438 static jobinfo_GeomOpt ji;
1439
1440 ji.prj = this;
1441 ji.go = geomopt_param(su); ji.go.Confirm();
1442 ji.show_dialog = false;
1443
1444 start_job_GeomOpt(& ji);
1445 return;
1446 }
1447
1448 if (!strcmp("mol_dyn", kw1)) // todo: how to set the options here?
1449 {
1450 setup * su = GetCurrentSetup();
1451 static jobinfo_MolDyn ji;
1452
1453 ji.prj = this;
1454 ji.md = moldyn_param(su); ji.md.Confirm();
1455 ji.show_dialog = false;
1456
1457 start_job_MolDyn(& ji);
1458 return;
1459 }
1460
1461 if (!strcmp("random_search", kw1))
1462 {
1463 char kw2[32]; istr >> kw2; // the 2nd keyword; cycles.
1464 char kw3[32]; istr >> kw3; // the 3rd keyword; optsteps.
1465
1466 static jobinfo_RandomSearch ji;
1467 char ** endptr = NULL;
1468
1469 ji.prj = this;
1470 ji.cycles = strtol(kw2, endptr, 10);
1471 ji.optsteps = strtol(kw3, endptr, 10);
1472
1473 start_job_RandomSearch(& ji);
1474 return;
1475 }
1476
1477 if (!strcmp("systematic_search", kw1))
1478 {
1479 char kw2[32]; istr >> kw2; // the 2nd keyword; divisions.
1480 char kw3[32]; istr >> kw3; // the 3rd keyword; optsteps.
1481
1482 char ** endptr = NULL;
1483 i32s divisions = strtol(kw2, endptr, 10);
1484 i32s optsteps = strtol(kw3, endptr, 10);
1485
1486 DoSystematicSearch(divisions, optsteps, true);
1487 return;
1488 }
1489
1490 if (!strcmp("montecarlo_search", kw1))
1491 {
1492 char kw2[32]; istr >> kw2; // the 2nd keyword; n_init_steps.
1493 char kw3[32]; istr >> kw3; // the 3rd keyword; n_simul_steps.
1494 char kw4[32]; istr >> kw4; // the 4th keyword; optsteps.
1495
1496 char ** endptr = NULL;
1497 i32s n_init_steps = strtol(kw2, endptr, 10);
1498 i32s n_simul_steps = strtol(kw3, endptr, 10);
1499 i32s optsteps = strtol(kw4, endptr, 10);
1500
1501 DoMonteCarloSearch(n_init_steps, n_simul_steps, optsteps, true);
1502 return;
1503 }
1504
1505 if (!strcmp("make_plot1", kw1))
1506 {
1507 char kw2[32]; istr >> kw2; // A
1508 char kw3[32]; istr >> kw3; // B
1509 char kw4[32]; istr >> kw4; // C
1510 char kw5[32]; istr >> kw5; // D
1511 char kw6[32]; istr >> kw6; // div
1512 char kw7[32]; istr >> kw7; // start_ang
1513 char kw8[32]; istr >> kw8; // end_ang
1514 char kw9[32]; istr >> kw9; // optsteps
1515 char ** endptr = NULL;
1516
1517 // show atom index 1,2,3,... to user ; it is 0,1,2,... internally!
1518 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1519
1520 i32s ia = strtol(kw2, endptr, 10) - 1;
1521 i32s ib = strtol(kw3, endptr, 10) - 1;
1522 i32s ic = strtol(kw4, endptr, 10) - 1;
1523 i32s id = strtol(kw5, endptr, 10) - 1;
1524 i32s div1 = strtol(kw6, endptr, 10);
1525 fGL start1 = strtod(kw7, endptr);
1526 fGL end1 = strtod(kw8, endptr);
1527
1528 i32s optsteps = strtol(kw9, endptr, 10);
1529
1530 DoEnergyPlot1D(ia, ib, ic, id, div1, start1, end1, optsteps);
1531 return;
1532 }
1533
1534 if (!strcmp("make_plot2", kw1))
1535 {
1536 char kw2[32]; istr >> kw2; // A
1537 char kw3[32]; istr >> kw3; // B
1538 char kw4[32]; istr >> kw4; // C
1539 char kw5[32]; istr >> kw5; // D
1540 char kw6[32]; istr >> kw6; // div
1541 char kw7[32]; istr >> kw7; // start_ang
1542 char kw8[32]; istr >> kw8; // end_ang
1543 char kw9[32]; istr >> kw9; // I
1544 char kwA[32]; istr >> kwA; // J
1545 char kwB[32]; istr >> kwB; // K
1546 char kwC[32]; istr >> kwC; // L
1547 char kwD[32]; istr >> kwD; // div
1548 char kwE[32]; istr >> kwE; // start_ang
1549 char kwF[32]; istr >> kwF; // end_ang
1550 char kwG[32]; istr >> kwG; // optsteps
1551 char ** endptr = NULL;
1552
1553 // show atom index 1,2,3,... to user ; it is 0,1,2,... internally!
1554 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1555
1556 i32s ia = strtol(kw2, endptr, 10) - 1;
1557 i32s ib = strtol(kw3, endptr, 10) - 1;
1558 i32s ic = strtol(kw4, endptr, 10) - 1;
1559 i32s id = strtol(kw5, endptr, 10) - 1;
1560 i32s div1 = strtol(kw6, endptr, 10);
1561 fGL start1 = strtod(kw7, endptr);
1562 fGL end1 = strtod(kw8, endptr);
1563
1564 i32s ii = strtol(kw9, endptr, 10) - 1;
1565 i32s ij = strtol(kwA, endptr, 10) - 1;
1566 i32s ik = strtol(kwB, endptr, 10) - 1;
1567 i32s il = strtol(kwC, endptr, 10) - 1;
1568 i32s div2 = strtol(kwD, endptr, 10);
1569 fGL start2 = strtod(kwE, endptr);
1570 fGL end2 = strtod(kwF, endptr);
1571
1572 i32s optsteps = strtol(kwG, endptr, 10);
1573
1574 DoEnergyPlot2D(ia, ib, ic, id, div1, start1, end1, ii, ij, ik, il, div2, start2, end2, optsteps);
1575 return;
1576 }
1577
1578 if (!strcmp("population_analysis_ESP", kw1))
1579 {
1580 setup1_qm * suqm = dynamic_cast<setup1_qm *>(current_setup);
1581 if (suqm == NULL) Message(_("Sorry, this is for QM models only!"));
1582 else
1583 {
1584 pop_ana_electrostatic pa(suqm);
1585 pa.DoPopAna();
1586
1587 // how to set the charge labels on in graphics?
1588 }
1589
1590 return;
1591 }
1592
1593 if (!strcmp("transition_state_search", kw1))
1594 {
1595 char kw2[32]; istr >> kw2; // the 2nd keyword; delta-E per step.
1596 char kw3[32]; istr >> kw3; // the 3rd keyword; initial force constant.
1597 char ** endptr = NULL;
1598
1599 fGL deltae = strtod(kw2, endptr);
1600 fGL initfc = strtod(kw3, endptr);
1601
1602 DoTransitionStateSearch(deltae, initfc);
1603 return;
1604 }
1605
1606 if (!strcmp("stationary_state_search", kw1))
1607 {
1608 char kw2[32]; istr >> kw2; // the 2nd keyword; steps.
1609 char ** endptr = NULL;
1610
1611 i32s steps = strtol(kw2, endptr, 10);
1612
1613 DoStationaryStateSearch(steps);
1614 return;
1615 }
1616
1617 if (!strcmp("set_current_orbital", kw1))
1618 {
1619 char kw2[32]; istr >> kw2; // the 2nd keyword; the orbital index.
1620
1621 char ** endptr = NULL;
1622 int index = strtol(kw2, endptr, 10);
1623 if (index < 0) index = 0;
1624
1625 qm_current_orbital = index;
1626
1627 ostringstream strR;
1628 strR << _("The current orbital is now ") << qm_current_orbital << "." << endl << ends;
1629 PrintToLog(strR.str().c_str());
1630 return;
1631 }
1632
1633 if (!strcmp("update_chains", kw1))
1634 {
1635 UpdateChains();
1636 PrintToLog("update_chains done.");
1637 return;
1638 }
1639
1640 if (!strcmp("build_amino", kw1))
1641 {
1642 char kw2[4096]; istr >> kw2; // sequence
1643 char kw3[32]; istr >> kw3; // helix/sheet (optional)
1644
1645 sb_chain_descriptor * sbcd = new sb_chain_descriptor(true);
1646
1647 for (int n1 = 0;n1 < strlen(kw2);n1++)
1648 {
1649 sbcd->AddRes1(kw2[n1]);
1650 }
1651
1652 // set only 3 torsions, and use defaults for the rest...
1653
1654 if (kw3[0] == 'h' || kw3[0] == 'H')
1655 {
1656 sbcd->def_tor.push_back(313.0 * M_PI / 180.0);
1657 sbcd->def_tor.push_back(180.0 * M_PI / 180.0);
1658 sbcd->def_tor.push_back(302.0 * M_PI / 180.0);
1659 }
1660 else
1661 {
1662 sbcd->def_tor.push_back(180.0 * M_PI / 180.0);
1663 sbcd->def_tor.push_back(180.0 * M_PI / 180.0);
1664 sbcd->def_tor.push_back(180.0 * M_PI / 180.0);
1665 }
1666
1667 model::amino_builder->Build(this, sbcd);
1668 delete sbcd;
1669
1670 UpdateAllGraphicsViews();
1671
1672 ostringstream strR;
1673 strR << _("built a sequence : ") << kw2 << endl << ends;
1674 PrintToLog(strR.str().c_str());
1675 return;
1676 }
1677
1678 if (!strcmp("build_nucleic", kw1))
1679 {
1680 char kw2[4096]; istr >> kw2; // sequence
1681
1682 sb_chain_descriptor * sbcd = new sb_chain_descriptor(true);
1683
1684 for (int n1 = 0;n1 < strlen(kw2);n1++)
1685 {
1686 sbcd->AddRes1(kw2[n1]);
1687 }
1688
1689 // set all 10 torsions...
1690
1691 sbcd->def_tor.push_back(261.0 * M_PI / 180.0);
1692 sbcd->def_tor.push_back(320.8 * M_PI / 180.0);
1693 sbcd->def_tor.push_back(208.6 * M_PI / 180.0);
1694 sbcd->def_tor.push_back(273.8 * M_PI / 180.0);
1695 sbcd->def_tor.push_back(105.6 * M_PI / 180.0);
1696 sbcd->def_tor.push_back(356.0 * M_PI / 180.0);
1697 sbcd->def_tor.push_back( 24.7 * M_PI / 180.0);
1698 sbcd->def_tor.push_back( 88.7 * M_PI / 180.0);
1699 sbcd->def_tor.push_back( 44.6 * M_PI / 180.0);
1700 sbcd->def_tor.push_back(264.6 * M_PI / 180.0);
1701
1702 model::nucleic_builder->Build(this, sbcd);
1703 delete sbcd;
1704
1705 UpdateAllGraphicsViews();
1706
1707 ostringstream strR;
1708 strR << _("built a sequence : ") << kw2 << endl << ends;
1709 PrintToLog(strR.str().c_str());
1710 return;
1711 }
1712
1713 if (!strcmp("orient", kw1))
1714 {
1715 fGL maxdim[3];
1716
1717 CenterCRDSet(0, true);
1718 OrientCRDSet(0, true, maxdim);
1719
1720 ostringstream strR;
1721 strR << _("maximum dimensions:") << endl << "X : " << maxdim[0] << endl;
1722 strR << "Y : " << maxdim[1] << endl << "Z : " << maxdim[2] << endl << ends;
1723 PrintToLog(strR.str().c_str());
1724
1725 UpdateAllGraphicsViews();
1726 return;
1727 }
1728
1729 if (!strcmp("solvate_box", kw1))
1730 {
1731 char kw2[32]; istr >> kw2; // xdim
1732 char kw3[32]; istr >> kw3; // ydim
1733 char kw4[32]; istr >> kw4; // zdim
1734 char kw5[32] = ""; istr >> kw5; // density (optional)
1735 char kw6[256] = ""; istr >> kw6; // filename (optional)
1736 char kw7[64] = ""; istr >> kw7; // export (optional)
1737
1738 char ** endptr = NULL;
1739 fGL xdim = strtod(kw2, endptr);
1740 fGL ydim = strtod(kw3, endptr);
1741 fGL zdim = strtod(kw4, endptr);
1742
1743 fGL density = 1.00; if (strlen(kw5) > 0) density = strtod(kw5, endptr);
1744 char * export_fn = NULL; if (!strcmp(kw7, "export")) export_fn = kw6;
1745
1746 dummy_project * solvent = NULL;
1747 if (strlen(kw6) > 0)
1748 {
1749 solvent = new dummy_project();
1750
1751 ostringstream fns;
1752 fns << kw6 << ".gpr" << ends;
1753
1754 ifstream ifile(fns.str().c_str(), ios::in);
1755 ReadGPR(* solvent, ifile, false);
1756 ifile.close();
1757 }
1758
1759 SolvateBox(xdim, ydim, zdim, density, solvent, export_fn);
1760 UpdateAllGraphicsViews();
1761 return;
1762 }
1763
1764 if (!strcmp("solvate_sphere", kw1))
1765 {
1766 char kw2[32]; istr >> kw2; // rad_solute
1767 char kw3[32]; istr >> kw3; // rad_solvent
1768 char kw4[32] = ""; istr >> kw4; // density (optional)
1769 char kw5[256] = ""; istr >> kw5; // filename (optional)
1770
1771 char ** endptr = NULL;
1772 // note 2005-01-02 : the MS compiler (vc7) choked on this
1773 // when the variable names rad1 and rad2 were used. HUH?!?!
1774 fGL radius1 = strtod(kw2, endptr);
1775 fGL radius2 = strtod(kw3, endptr);
1776
1777 fGL density = 1.00; // in kg/liter as usual...
1778 if (strlen(kw4) > 0) density = strtod(kw4, endptr);
1779
1780 dummy_project * solvent = NULL;
1781 if (strlen(kw5) > 0)
1782 {
1783 solvent = new dummy_project();
1784
1785 ostringstream fns;
1786 fns << kw5 << ".gpr" << ends;
1787
1788 ifstream ifile(fns.str().c_str(), ios::in);
1789 ReadGPR(* solvent, ifile, false);
1790 ifile.close();
1791 }
1792
1793 SolvateSphere(radius1, radius2, density, solvent);
1794 UpdateAllGraphicsViews();
1795 return;
1796 }
1797
1798 if (!strcmp("set_formal_charge", kw1))
1799 {
1800 char kw2[32]; istr >> kw2; // the 2nd keyword; index.
1801 char kw3[32]; istr >> kw3; // the 3rd keyword; charge.
1802 char ** endptr = NULL;
1803
1804 i32s index = strtol(kw2, endptr, 10);
1805 i32s charge = strtol(kw3, endptr, 10);
1806
1807 if (!IsIndexClean()) UpdateIndex();
1808
1809 atom * atmr = NULL;
1810 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
1811 {
1812 if ((* it1).index == index) { atmr = & (* it1); break; }
1813 }
1814
1815 if (atmr != NULL)
1816 {
1817 atmr->formal_charge = charge;
1818 UpdateAllGraphicsViews(); // update the labels...
1819 }
1820 else
1821 {
1822 ostringstream strE;
1823 strE << _("Sorry, atom not found!") << endl << ends;
1824 PrintToLog(strE.str().c_str());
1825 }
1826 return;
1827 }
1828
1829 if (!strcmp("evaluate_Bfact", kw1))
1830 {
1831 EvaluateBFact();
1832 return;
1833 }
1834
1835 if (!strcmp("evaluate_diffconst", kw1))
1836 {
1837 char kw2[32]; istr >> kw2; // the 2nd keyword; dt.
1838
1839 char ** endptr = NULL;
1840 fGL dt = strtod(kw2, endptr);
1841
1842 EvaluateDiffConst(dt);
1843 return;
1844 }
1845
1846 // if the command is not recognized above, we will print out an error message here.
1847
1848 ostringstream strE;
1849 strE << _("ERROR : Unknown command : ") << command << endl;
1850 strE << _("The \"help\" command will give more information about command strings.") << endl;
1851 strE << ends;
1852
1853 PrintToLog(strE.str().c_str());
1854 }
1855
1856 /*##############################################*/
1857 /*##############################################*/
1858
DoDeleteCurrentObject(void)1859 void project::DoDeleteCurrentObject(void)
1860 {
1861 if (selected_object != NULL)
1862 {
1863 bool test1 = base_app::GetAppB()->RemoveLight(selected_object);
1864 bool test2 = test1; if (!test1) test2 = RemoveObject(selected_object);
1865
1866 if (test2)
1867 {
1868 selected_object = NULL;
1869 UpdateAllGraphicsViews();
1870 }
1871 }
1872 }
1873
DoSwitchLocalLights(ogl_camera * cam,bool report)1874 void project::DoSwitchLocalLights(ogl_camera * cam, bool report)
1875 {
1876 cam->use_local_lights = !cam->use_local_lights;
1877 if (report) cout << _("local lights = ") << (cam->use_local_lights ? _("on") : _("off")) << endl;
1878 base_app::GetAppB()->SetupLights(cam); UpdateGraphicsViews(cam);
1879 }
1880
DoSwitchGlobalLights(ogl_camera * cam,bool report)1881 void project::DoSwitchGlobalLights(ogl_camera * cam, bool report)
1882 {
1883 cam->use_global_lights = !cam->use_global_lights;
1884 if (report) cout << _("global lights = ") << (cam->use_global_lights ? _("on") : _("off")) << endl;
1885 base_app::GetAppB()->SetupLights(cam); UpdateGraphicsViews(cam);
1886 }
1887
GetDefaultFocus(void)1888 fGL project::GetDefaultFocus(void)
1889 {
1890 return 2.0;
1891 }
1892
GetDefaultColorMode(void)1893 color_mode * project::GetDefaultColorMode(void)
1894 {
1895 return & project::cm_element;
1896 }
1897
SelectAll(void)1898 void project::SelectAll(void)
1899 {
1900 if (selected_object != NULL)
1901 {
1902 selected_object = NULL;
1903 }
1904
1905 iter_al it1 = atom_list.begin();
1906 while (it1 != atom_list.end()) (* it1++).flags |= ATOMFLAG_USER_SELECTED;
1907
1908 UpdateAllGraphicsViews();
1909 }
1910
InvertSelection(void)1911 void project::InvertSelection(void)
1912 {
1913 if (selected_object != NULL)
1914 {
1915 selected_object = NULL;
1916 }
1917
1918 iter_al it1 = atom_list.begin();
1919 while (it1 != atom_list.end()) (* it1++).flags ^= ATOMFLAG_USER_SELECTED;
1920
1921 UpdateAllGraphicsViews();
1922 }
1923
HideSelected(void)1924 void project::HideSelected(void)
1925 {
1926 iter_al it1 = atom_list.begin();
1927 while (it1 != atom_list.end())
1928 {
1929 if ((* it1).flags & ATOMFLAG_USER_SELECTED)
1930 {
1931 (* it1).flags |= ATOMFLAG_USER_HIDDEN;
1932 }
1933
1934 it1++;
1935 }
1936
1937 UpdateAllGraphicsViews();
1938 }
1939
ShowSelected(void)1940 void project::ShowSelected(void)
1941 {
1942 iter_al it1 = atom_list.begin();
1943 while (it1 != atom_list.end())
1944 {
1945 if ((* it1).flags & ATOMFLAG_USER_SELECTED)
1946 {
1947 (* it1).flags &= (~ATOMFLAG_USER_HIDDEN);
1948 }
1949
1950 it1++;
1951 }
1952
1953 UpdateAllGraphicsViews();
1954 }
1955
LockSelected(void)1956 void project::LockSelected(void)
1957 {
1958 iter_al it1 = atom_list.begin();
1959 while (it1 != atom_list.end())
1960 {
1961 if ((* it1).flags & ATOMFLAG_USER_SELECTED)
1962 {
1963 (* it1).flags |= ATOMFLAG_USER_LOCKED;
1964 }
1965
1966 it1++;
1967 }
1968
1969 UpdateAllGraphicsViews();
1970 }
1971
UnlockSelected(void)1972 void project::UnlockSelected(void)
1973 {
1974 iter_al it1 = atom_list.begin();
1975 while (it1 != atom_list.end())
1976 {
1977 if ((* it1).flags & ATOMFLAG_USER_SELECTED)
1978 {
1979 (* it1).flags &= (~ATOMFLAG_USER_LOCKED);
1980 }
1981
1982 it1++;
1983 }
1984
1985 UpdateAllGraphicsViews();
1986 }
1987
DeleteSelected(void)1988 void project::DeleteSelected(void)
1989 {
1990 if (selected_object != NULL)
1991 {
1992 selected_object = NULL; // is this right??? memoryleak at least...
1993 }
1994
1995 iter_al it1 = atom_list.begin();
1996 while (it1 != atom_list.end())
1997 {
1998 if ((* it1).flags & ATOMFLAG_USER_SELECTED)
1999 {
2000 RemoveAtom(it1);
2001 it1 = atom_list.begin(); // reset the search!!!
2002 }
2003 else it1++;
2004 }
2005
2006 UpdateAllGraphicsViews();
2007 }
2008
TestAtom(atom * ref,rmode rm)2009 bool project::TestAtom(atom * ref, rmode rm)
2010 {
2011 if (ref->flags & (ATOMFLAG_IS_HIDDEN | ATOMFLAG_USER_HIDDEN)) return false;
2012
2013 if (rm == Transform1 && (ref->flags & ATOMFLAG_USER_SELECTED)) return false;
2014 if (rm == Transform2 && !(ref->flags & ATOMFLAG_USER_SELECTED)) return false;
2015
2016 return true;
2017 }
2018
TestBond(bond * ref,rmode rm)2019 bool project::TestBond(bond * ref, rmode rm)
2020 {
2021 if (ref->atmr[0]->flags & (ATOMFLAG_IS_HIDDEN | ATOMFLAG_USER_HIDDEN)) return false;
2022 if (ref->atmr[1]->flags & (ATOMFLAG_IS_HIDDEN | ATOMFLAG_USER_HIDDEN)) return false;
2023
2024 if (rm == Transform1 && (ref->atmr[0]->flags & ATOMFLAG_USER_SELECTED)) return false; // no need to study the another...
2025 if (rm == Transform2 && !(ref->atmr[0]->flags & ATOMFLAG_USER_SELECTED)) return false; // ...atom due to the test below?
2026
2027 bool test1 = (ref->atmr[0]->flags & ATOMFLAG_USER_SELECTED) ? true : false;
2028 bool test2 = (ref->atmr[1]->flags & ATOMFLAG_USER_SELECTED) ? true : false;
2029 if (rm != Normal && test1 != test2) return false;
2030
2031 if (ref->do_not_render_TSS_fixmelater) return false; // temporary, for transition_state_search only...
2032
2033 return true;
2034 }
2035
SetColor(color_mode * cm,atom * ref,bool black_and_white)2036 void project::SetColor(color_mode * cm, atom * ref, bool black_and_white)
2037 {
2038 fGL select_color[3] = { 1.0, 0.0, 1.0 };
2039 fGL measure_color[3] = { 0.0, 1.0, 1.0 };
2040
2041 if (ref->flags & (ATOMFLAG_USER_SELECTED | ATOMFLAG_MEASURE_TOOL_SEL))
2042 {
2043 fGL * color = select_color;
2044 if (ref->flags & ATOMFLAG_MEASURE_TOOL_SEL) color = measure_color;
2045
2046 if (black_and_white) // if we have a red/blue stereo mode, average the colours to shades of gray!
2047 {
2048 fGL average = (color[0] + color[1] + color[2]) / 3.0;
2049 color[0] = color[1] = color[2] = average;
2050 }
2051
2052 glColor3f(color[0], color[1], color[2]);
2053 }
2054 else
2055 {
2056 static fGL color[4];
2057 cm->GetColor4(ref, -1, color);
2058
2059 if (black_and_white) // if we have a red/blue stereo mode, average the colours to shades of gray!
2060 {
2061 fGL average = (color[0] + color[1] + color[2]) / 3.0;
2062 color[0] = color[1] = color[2] = average;
2063 }
2064
2065 glColor3fv(color);
2066 }
2067 }
2068
DrawCylinder1(const fGL ** crd,const fGL ** col,const fGL * rad)2069 void project::DrawCylinder1(const fGL ** crd, const fGL ** col, const fGL * rad)
2070 {
2071 fGL rsum = rad[0] + rad[1];
2072
2073 for (i32s n1 = 0;n1 < 2;n1++)
2074 {
2075 glColor3fv(col[n1]);
2076
2077 v3d<fGL> crt = v3d<fGL>(crd[n1], crd[!n1]);
2078 fGL pol[3]; crt2pol(crt.data, pol);
2079
2080 const int resolution = 10;
2081
2082 GLUquadricObj * qo = gluNewQuadric();
2083 gluQuadricDrawStyle(qo, (GLenum) GLU_FILL); glPushMatrix();
2084
2085 glTranslated(crd[n1][0], crd[n1][1], crd[n1][2]);
2086
2087 glRotated(180.0 * pol[1] / M_PI, 0.0, 1.0, 0.0);
2088 glRotated(180.0 * pol[2] / M_PI, sin(-pol[1]), 0.0, cos(-pol[1]));
2089
2090 fGL length = crt.len() * rad[n1] / rsum;
2091 gluCylinder(qo, 0.1*rad[n1], 0.1*rad[!n1], length, resolution, resolution / 2);
2092
2093 glPopMatrix(); gluDeleteQuadric(qo);
2094 }
2095 }
2096
Render(oglview_wcl * wcl,rmode rm)2097 void project::Render(oglview_wcl * wcl, rmode rm)
2098 {
2099 const fGL label_color[3] = { 0.0, 1.0, 1.0 }; // looks bad but won't fade easily into other colours...
2100
2101 bool accum = wcl->accumulate; if (rm != Normal) accum = false;
2102 //if (accum) { glClear(GL_ACCUM_BUFFER_BIT); UpdateAccumValues(); }
2103 //else if (rm != Transform2) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
2104
2105 if (use_boundary_potential && rm == Normal)
2106 {
2107 for (int loop = 0;loop < 2;loop++)
2108 {
2109 fGL radius;
2110
2111 if (!loop)
2112 {
2113 glColor3f(0.5, 0.0, 0.5);
2114 radius = saved_boundary_potential_rad_solute;
2115 }
2116 else
2117 {
2118 glColor3f(1.0, 0.0, 1.0);
2119 radius = saved_boundary_potential_rad_solvent;
2120 }
2121
2122 glPushMatrix();
2123 glTranslated(0.0, 0.0, 0.0); // TODO : set the engine::bp_center[] here!!!
2124
2125 glBegin(GL_LINES);
2126
2127 fGL ang1 = 0.0;
2128 fGL ca1 = radius * cos(ang1);
2129 fGL sa1 = radius * sin(ang1);
2130
2131 const i32s divisions = 12;
2132 for (i32s n1 = 0;n1 < divisions;n1++)
2133 {
2134 fGL ang2 = 2.0 * M_PI * ((fGL) (n1 + 1)) / (fGL) divisions;
2135 fGL ca2 = radius * cos(ang2);
2136 fGL sa2 = radius * sin(ang2);
2137
2138 glVertex3f(ca1, sa1, 0.0);
2139 glVertex3f(ca2, sa2, 0.0);
2140
2141 glVertex3f(ca1, 0.0, sa1);
2142 glVertex3f(ca2, 0.0, sa2);
2143
2144 glVertex3f(0.0, ca1, sa1);
2145 glVertex3f(0.0, ca2, sa2);
2146
2147 ang1 = ang2; ca1 = ca2; sa1 = sa2;
2148 }
2149
2150 glEnd();
2151
2152 glPopMatrix();
2153 }
2154 }
2155
2156 if (use_periodic_boundary_conditions && rm == Normal)
2157 {
2158 glLineWidth(1.0);
2159 glColor3f(1.0, 0.0, 1.0);
2160 glBegin(GL_LINES);
2161
2162 glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2163 glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2164
2165 glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2166 glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2167
2168 glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2169 glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2170
2171 glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2172 glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2173
2174 glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2175 glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2176
2177 glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2178 glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2179
2180 glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2181 glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2182
2183 glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2184 glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2185
2186 glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2187 glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2188
2189 glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2190 glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2191
2192 glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2193 glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
2194
2195 glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2196 glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
2197
2198 glEnd();
2199 }
2200
2201 if (wcl->enable_fog) glEnable(GL_FOG);
2202
2203 i32s layers = 0;
2204 //if (use_periodic_boundary_conditions && rm == Normal) layers = 1; // un-comment this to render the periodic images...
2205
2206 for (i32s r1 = -layers;r1 < (layers + 1);r1++)
2207 {
2208 for (i32s r2 = -layers;r2 < (layers + 1);r2++)
2209 {
2210 for (i32s r3 = -layers;r3 < (layers + 1);r3++)
2211 {
2212 glPushMatrix();
2213
2214 fGL trans1 = r1 * 2.0 * saved_periodic_box_HALFdim[0];
2215 fGL trans2 = r2 * 2.0 * saved_periodic_box_HALFdim[1];
2216 fGL trans3 = r3 * 2.0 * saved_periodic_box_HALFdim[2];
2217
2218 glTranslated(trans1, trans2, trans3);
2219
2220 RenderOnce(wcl, rm, accum);
2221
2222 glPopMatrix();
2223 }
2224 }
2225 }
2226
2227 if (accum) glAccum(GL_RETURN, 1.0);
2228 else if (rm != Transform2) RenderObjects(wcl);
2229
2230 if (wcl->label == LABEL_INDEX)
2231 {
2232 // show atom index 1,2,3,... to user ; it is 0,1,2,... internally!
2233 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2234
2235 i32s tmp1 = 1;
2236
2237 glColor3f(label_color[0], label_color[1], label_color[2]);
2238 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
2239 {
2240 if ((* it1).flags & ATOMFLAG_IS_HIDDEN) { tmp1++; continue; }
2241
2242 ostringstream str;
2243 str << tmp1++ << ends;
2244
2245 const fGL * cdata = (* it1).GetCRD(0);
2246 fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];
2247
2248 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2249 }
2250 }
2251 else if (wcl->label == LABEL_F_CHARGE)
2252 {
2253 glColor3f(label_color[0], label_color[1], label_color[2]);
2254 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
2255 {
2256 if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
2257
2258 ostringstream str;
2259 str.setf(ios::fixed | ios::showpos); str << (* it1).formal_charge << ends;
2260
2261 const fGL * cdata = (* it1).GetCRD(0);
2262 fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];
2263
2264 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2265 }
2266 }
2267 else if (wcl->label == LABEL_P_CHARGE)
2268 {
2269 glColor3f(label_color[0], label_color[1], label_color[2]);
2270 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
2271 {
2272 if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
2273
2274 ostringstream str;
2275 str.setf(ios::fixed | ios::showpos); str.precision(4); str << (* it1).charge << ends;
2276
2277 const fGL * cdata = (* it1).GetCRD(0);
2278 fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];
2279
2280 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2281 }
2282 }
2283 else if (wcl->label == LABEL_ELEMENT)
2284 {
2285 glColor3f(label_color[0], label_color[1], label_color[2]);
2286 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
2287 {
2288 if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
2289
2290 ostringstream str;
2291 str << (* it1).el.GetSymbol() << ends;
2292
2293 const fGL * cdata = (* it1).GetCRD(0);
2294 fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];
2295
2296 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2297 }
2298 }
2299 else if (wcl->label == LABEL_BUILDER_ID)
2300 {
2301 glColor3f(label_color[0], label_color[1], label_color[2]);
2302 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
2303 {
2304 if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
2305
2306 ostringstream str;
2307 str << "0x" << hex << (* it1).builder_res_id << ends;
2308
2309 const fGL * cdata = (* it1).GetCRD(0);
2310 fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];
2311
2312 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2313 }
2314 }
2315 else if (wcl->label == LABEL_ATOMTYPE)
2316 {
2317 glColor3f(label_color[0], label_color[1], label_color[2]);
2318 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
2319 {
2320 if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
2321
2322 ostringstream str;
2323 str << "0x" << hex << (* it1).atmtp << ends;
2324
2325 // str << "0x" << hex << (* it1).atmtp_E << ends; // debug...
2326 // str << (* it1).atmtp_s << ends; // debug...
2327
2328 /* if (!(* it1).atRS) str << "none" << ends;
2329 else
2330 {
2331 atomtype_mmRS * atmtp = (* it1).atRS;
2332 if (!atmtp) str << "0x" << hex << (* it1).atmtp << ends;
2333
2334 for (int n1 = 0;n1 < atmtp->GetSize();n1++)
2335 {
2336 str << "0x" << hex << atmtp->GetAtomType(n1) << " (" << atmtp->GetWeight(n1) << ")" << endl;
2337 } str << ends;
2338 } */
2339
2340 const fGL * cdata = (* it1).GetCRD(0);
2341 fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];
2342
2343 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2344 }
2345 }
2346 else if (wcl->label == LABEL_BONDTYPE)
2347 {
2348 glColor3f(label_color[0], label_color[1], label_color[2]);
2349 for (iter_bl it1 = bond_list.begin();it1 != bond_list.end();it1++)
2350 {
2351 if ((* it1).atmr[0]->flags & ATOMFLAG_IS_HIDDEN) continue;
2352 if ((* it1).atmr[1]->flags & ATOMFLAG_IS_HIDDEN) continue;
2353
2354 ostringstream str;
2355 str << (* it1).bt.GetSymbol1() << ends;
2356
2357 const fGL * cd1 = (* it1).atmr[0]->GetCRD(0); const fGL * cd2 = (* it1).atmr[1]->GetCRD(0);
2358 fGL x = (cd1[0] + cd2[0]) / 2.0; fGL y = (cd1[1] + cd2[1]) / 2.0; fGL z = (cd1[2] + cd2[2]) / 2.0;
2359
2360 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2361 }
2362 }
2363 else if (wcl->label == LABEL_RESIDUE)
2364 {
2365 if (ref_civ != NULL)
2366 {
2367 glColor3f(label_color[0], label_color[1], label_color[2]);
2368
2369 vector<chn_info> & ci_vector = (* ref_civ);
2370 for (i32u chn = 0;chn < ci_vector.size();chn++)
2371 {
2372 iter_al range1[2]; GetRange(1, chn, range1);
2373 const char * tmp_seq1 = ci_vector[chn].GetSequence1();
2374
2375 for (i32s res = 0;res < ci_vector[chn].GetLength();res++)
2376 {
2377 iter_al range2[2]; GetRange(2, range1, res, range2);
2378 fGL rescrd[3] = { 0.0, 0.0, 0.0 }; i32s counter = 0;
2379
2380 // SLOW because coordinates calculated on-the-fly!!! save them somewhere???
2381 // SLOW because coordinates calculated on-the-fly!!! save them somewhere???
2382 // SLOW because coordinates calculated on-the-fly!!! save them somewhere???
2383
2384 for (iter_al it1 = range2[0];it1 != range2[1];it1++)
2385 {
2386 const fGL * atmcrd = (* it1).GetCRD(0);
2387 rescrd[0] += atmcrd[0]; rescrd[1] += atmcrd[1]; rescrd[2] += atmcrd[2];
2388 counter++;
2389 }
2390
2391 fGL x = rescrd[0] / (fGL) counter;
2392 fGL y = rescrd[1] / (fGL) counter;
2393 fGL z = rescrd[2] / (fGL) counter;
2394
2395 // show chn/res index 1,2,3,... to user ; it is 0,1,2,... internally!
2396 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2397
2398 ostringstream str;
2399 str << tmp_seq1[res] << " (" << (chn + 1) << "/" << (res + 1) << ")" << ends;
2400
2401 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2402 }
2403 }
2404 }
2405 }
2406 else if (wcl->label == LABEL_SEC_STRUCT)
2407 {
2408 if (ref_civ != NULL)
2409 {
2410 glColor3f(label_color[0], label_color[1], label_color[2]);
2411
2412 vector<chn_info> & ci_vector = (* ref_civ);
2413 for (i32u chn = 0;chn < ci_vector.size();chn++)
2414 {
2415 iter_al range1[2]; GetRange(1, chn, range1);
2416 const char * tmp_states = ci_vector[chn].GetSecStrStates();
2417
2418 for (i32s res = 0;res < ci_vector[chn].GetLength();res++)
2419 {
2420 iter_al range2[2]; GetRange(2, range1, res, range2);
2421 fGL rescrd[3] = { 0.0, 0.0, 0.0 }; i32s counter = 0;
2422
2423 // SLOW because coordinates calculated on-the-fly!!! save them somewhere???
2424 // SLOW because coordinates calculated on-the-fly!!! save them somewhere???
2425 // SLOW because coordinates calculated on-the-fly!!! save them somewhere???
2426
2427 for (iter_al it1 = range2[0];it1 != range2[1];it1++)
2428 {
2429 const fGL * atmcrd = (* it1).GetCRD(0);
2430 rescrd[0] += atmcrd[0]; rescrd[1] += atmcrd[1]; rescrd[2] += atmcrd[2];
2431 counter++;
2432 }
2433
2434 fGL x = rescrd[0] / (fGL) counter;
2435 fGL y = rescrd[1] / (fGL) counter;
2436 fGL z = rescrd[2] / (fGL) counter;
2437
2438 ostringstream str;
2439 str << tmp_states[res] << ends;
2440
2441 wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
2442 }
2443 }
2444 }
2445 }
2446
2447 if (wcl->enable_fog) glDisable(GL_FOG);
2448
2449 // finally call this to handle transparency...
2450 // finally call this to handle transparency...
2451 // finally call this to handle transparency...
2452
2453 base_app::GetAppB()->RenderAllTPs(wcl->GetCam());
2454 }
2455
RenderOnce(oglview_wcl * wcl,rmode rm,bool accum)2456 void project::RenderOnce(oglview_wcl * wcl, rmode rm, bool accum)
2457 {
2458 bool do_bw = (wcl->GetCam()->stereo_mode && !wcl->GetCam()->stereo_relaxed);
2459
2460 for (i32u n1 = 0;n1 < cs_vector.size();n1++)
2461 {
2462 if (!GetCRDSetVisible(n1)) continue;
2463 if (accum) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); // FIXME!!!
2464
2465 if (wcl->render == RENDER_WIREFRAME)
2466 {
2467 glPointSize(3.0); glLineWidth(1.0);
2468 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++) // wireframe atoms
2469 {
2470 if (!TestAtom(& (* it1), rm)) continue;
2471 glPushName(GLNAME_MD_TYPE1); glPushName((* it1).my_glname);
2472
2473 glBegin(GL_POINTS);
2474 SetColor(wcl->colormode, & (* it1), do_bw);
2475 glVertex3fv((* it1).GetCRD(n1));
2476 glEnd();
2477
2478 glPopName(); glPopName();
2479 }
2480
2481 glEnable(GL_LINE_STIPPLE);
2482 for (iter_bl it2 = bond_list.begin();it2 != bond_list.end();it2++) // wireframe bonds
2483 {
2484 if (!TestBond(& (* it2), rm)) continue;
2485
2486 switch ((* it2).bt.GetSymbol1())
2487 {
2488 case 'S': glLineStipple(1, 0xFFFF); break;
2489 case 'C': glLineStipple(1, 0x3FFF); break;
2490 case 'D': glLineStipple(1, 0x3F3F); break;
2491 case 'T': glLineStipple(1, 0x3333); break;
2492 }
2493
2494 glBegin(GL_LINES);
2495 SetColor(wcl->colormode, (* it2).atmr[0], do_bw);
2496 glVertex3fv((* it2).atmr[0]->GetCRD(n1));
2497 SetColor(wcl->colormode, (* it2).atmr[1], do_bw);
2498 glVertex3fv((* it2).atmr[1]->GetCRD(n1));
2499 glEnd();
2500 }
2501 glDisable(GL_LINE_STIPPLE);
2502 }
2503
2504 if (wcl->render != RENDER_WIREFRAME && wcl->render != RENDER_NOTHING)
2505 {
2506 glEnable(GL_LIGHTING);
2507
2508 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++) // atoms as spheres
2509 {
2510 if (!TestAtom(& (* it1), rm)) continue;
2511
2512 SetColor(wcl->colormode, & (* it1), do_bw);
2513
2514 float rad = 0.0; int res = 0;
2515 switch (wcl->render)
2516 {
2517 case RENDER_BALL_AND_STICK:
2518 rad = 0.035;
2519 res = 12;
2520 break;
2521
2522 case RENDER_VAN_DER_WAALS:
2523 rad = (* it1).vdwr;
2524 res = 22;
2525 break;
2526
2527 case RENDER_CYLINDERS:
2528 rad = 0.035;
2529 res = 12;
2530 break;
2531 }
2532
2533 glPushName(GLNAME_MD_TYPE1); glPushName((* it1).my_glname);
2534
2535 GLUquadricObj * qo = gluNewQuadric();
2536 gluQuadricDrawStyle(qo, (GLenum) GLU_FILL);
2537
2538 glPushMatrix();
2539 const fGL * cdata = (* it1).GetCRD(n1);
2540 glTranslated(cdata[0], cdata[1], cdata[2]);
2541 gluSphere(qo, rad, res, res / 2);
2542 glPopMatrix();
2543 gluDeleteQuadric(qo);
2544
2545 glPopName(); glPopName();
2546 }
2547
2548 glDisable(GL_LIGHTING);
2549 }
2550
2551 if (wcl->render == RENDER_BALL_AND_STICK || wcl->render == RENDER_CYLINDERS)
2552 {
2553 glEnable(GL_LIGHTING);
2554
2555 for (iter_bl it1 = bond_list.begin();it1 != bond_list.end();it1++) // bonds as cylinders
2556 {
2557 if (!TestBond(& (* it1), rm)) continue;
2558
2559 fGL vdwr[2] =
2560 {
2561 (* it1).atmr[0]->vdwr,
2562 (* it1).atmr[1]->vdwr
2563 };
2564
2565 fGL vdwrsum = vdwr[0] + vdwr[1];
2566
2567 for (i32s n2 = 0;n2 < 2;n2++)
2568 {
2569 const fGL * crd1 = (* it1).atmr[n2]->GetCRD(n1);
2570 const fGL * crd2 = (* it1).atmr[!n2]->GetCRD(n1);
2571 v3d<fGL> crt1 = v3d<fGL>(crd1);
2572 v3d<fGL> crt2 = v3d<fGL>(crd2);
2573 v3d<fGL> crt = crt2 - crt1;
2574
2575 fGL pol[3]; crt2pol(crt.data, pol);
2576
2577 SetColor(wcl->colormode, (* it1).atmr[n2], do_bw);
2578
2579 float trans, rad = 0.0; int res = 0;
2580 switch (wcl->render)
2581 {
2582 case RENDER_BALL_AND_STICK:
2583 rad = 0.01;
2584 res = 6;
2585 break;
2586
2587 case RENDER_CYLINDERS:
2588 rad = 0.035;
2589 res = 12;
2590 break;
2591 }
2592
2593 glPushName(GLNAME_MD_TYPE1); glPushName((* it1).atmr[n2]->my_glname);
2594
2595 GLUquadricObj * qo = gluNewQuadric();
2596 gluQuadricDrawStyle(qo, (GLenum) GLU_FILL);
2597 glPushMatrix();
2598
2599 glTranslated(crd1[0], crd1[1], crd1[2]);
2600
2601 glRotated(180.0 * pol[1] / M_PI, 0.0, 1.0, 0.0);
2602 glRotated(180.0 * pol[2] / M_PI, sin(-pol[1]), 0.0, cos(-pol[1]));
2603
2604 // any chance to further define the orientation of, for example, double bonds???
2605 // one more rotation would be needed. but what is the axis, and how much to rotate???
2606
2607 fGL length = crt.len() * vdwr[n2] / vdwrsum;
2608
2609 if (wcl->render == RENDER_BALL_AND_STICK)
2610 switch ((* it1).bt.GetValue())
2611 {
2612 case BONDTYPE_DOUBLE:
2613 trans = rad;
2614 rad = rad / 1.5;
2615
2616 if (n2)
2617 glTranslated(0.0, trans, 0.0);
2618 else
2619 glTranslated(0.0, -trans, 0.0);
2620 gluCylinder(qo, rad, rad, length, res, 1);
2621 if (n2)
2622 glTranslated(0.0, -2.0 * trans, 0.0);
2623 else
2624 glTranslated(0.0, 2.0 * trans, 0.0);
2625 gluCylinder(qo, rad, rad, length, res, 1);
2626 break;
2627
2628 case BONDTYPE_CNJGTD:
2629 trans = rad;
2630 rad = rad / 1.5;
2631
2632 if (n2)
2633 glTranslated(0.0, trans, 0.0);
2634 else
2635 glTranslated(0.0, -trans, 0.0);
2636 gluCylinder(qo, rad, rad, length, res, 1);
2637 if (n2)
2638 glTranslated(0.0, -2.0 * trans, 0.0);
2639 else
2640 glTranslated(0.0, 2.0 * trans, 0.0);
2641
2642 glEnable(GL_LINE_STIPPLE);
2643 glLineStipple(1, 0x3F3F);
2644 gluQuadricDrawStyle(qo, (GLenum) GLU_LINE);
2645 gluCylinder(qo, rad, rad, length, res, 1);
2646 glDisable(GL_LINE_STIPPLE);
2647 break;
2648
2649 case BONDTYPE_TRIPLE:
2650 trans = rad;
2651 rad = rad / 2.0;
2652
2653 if (n2)
2654 glTranslated(0.0, trans, 0.0);
2655 else
2656 glTranslated(0.0, -trans, 0.0);
2657 gluCylinder(qo, rad, rad, length, res, 1);
2658 if (n2)
2659 glTranslated(0.0, -trans, 0.0);
2660 else
2661 glTranslated(0.0, trans, 0.0);
2662 gluCylinder(qo, rad, rad, length, res, 1);
2663 if (n2)
2664 glTranslated(0.0, -trans, 0.0);
2665 else
2666 glTranslated(0.0, trans, 0.0);
2667 gluCylinder(qo, rad, rad, length, res,1);
2668 break;
2669
2670 default:
2671 gluCylinder(qo, rad, rad, length, res, 1);
2672 }
2673 else
2674 gluCylinder(qo, rad, rad, length, res, 1);
2675
2676 glPopMatrix();
2677 gluDeleteQuadric(qo);
2678
2679 glPopName(); glPopName();
2680 }
2681 }
2682
2683 glDisable(GL_LIGHTING);
2684 }
2685
2686 // render the additional stuff related to SF.
2687 // render the additional stuff related to SF.
2688 // render the additional stuff related to SF.
2689
2690 setup1_sf * susf = dynamic_cast<setup1_sf *>(current_setup);
2691 if (susf != NULL)
2692 {
2693 if (susf->mode == setup1_sf::modeP3)
2694 {
2695 for (i32u n2 = 0;n2 < susf->hi_vector.size();n2++) // helix constraints
2696 {
2697 glEnable(GL_LINE_STIPPLE);
2698 glLineStipple(1, 0x3333);
2699
2700 glBegin(GL_LINES);
2701 glColor3f(0.20, 1.00, 0.10);
2702
2703 const int sz = susf->hi_vector[n2].ca_H_don.size();
2704 if (sz != (int) susf->hi_vector[n2].ca_H_acc.size()) assertion_failed(__FILE__, __LINE__, "donHmc/accHmc mismatch");
2705
2706 for (i32u n3 = 0;n3 < sz;n3++)
2707 {
2708 const fGL * crd1 = susf->hi_vector[n2].ca_H_don[n3]->GetCRD(n1);
2709 const fGL * crd2 = susf->hi_vector[n2].ca_H_acc[n3]->GetCRD(n1);
2710
2711 glVertex3fv(crd1);
2712 glVertex3fv(crd2);
2713 }
2714
2715 glEnd();
2716 glDisable(GL_LINE_STIPPLE);
2717 }
2718
2719 for (i32u n2 = 0;n2 < susf->sp_vector.size();n2++) // strand constraints ; straight
2720 {
2721 glEnable(GL_LINE_STIPPLE);
2722 glLineStipple(1, 0x3333);
2723
2724 glBegin(GL_LINES);
2725 glColor3f(0.20, 1.00, 0.10);
2726
2727 const int sz = susf->sp_vector[n2].ca_S_2x.size();
2728 if (sz % 2 != 0) assertion_failed(__FILE__, __LINE__, "ca_S_2x has an odd size.");
2729
2730 for (i32u n3 = 0;n3 < sz;n3 += 2)
2731 {
2732 const fGL * crd1 = susf->sp_vector[n2].ca_S_2x[n3 + 0]->GetCRD(n1);
2733 const fGL * crd2 = susf->sp_vector[n2].ca_S_2x[n3 + 1]->GetCRD(n1);
2734
2735 glVertex3fv(crd1);
2736 glVertex3fv(crd2);
2737 }
2738
2739 glEnd();
2740 glDisable(GL_LINE_STIPPLE);
2741 }
2742
2743 for (i32u n2 = 0;n2 < susf->sp_vector.size();n2++) // strand constraints ; crossed
2744 {
2745 glEnable(GL_LINE_STIPPLE);
2746 glLineStipple(1, 0x3333);
2747
2748 glBegin(GL_LINES);
2749 glColor3f(0.80, 1.00, 0.10);
2750
2751 const int sz = susf->sp_vector[n2].cx_S_2x.size();
2752 if (sz % 2 != 0) assertion_failed(__FILE__, __LINE__, "cx_S_2x has an odd size.");
2753
2754 for (i32u n3 = 0;n3 < sz;n3 += 2)
2755 {
2756 const fGL * crd1 = susf->sp_vector[n2].cx_S_2x[n3 + 0]->GetCRD(n1);
2757 const fGL * crd2 = susf->sp_vector[n2].cx_S_2x[n3 + 1]->GetCRD(n1);
2758
2759 glVertex3fv(crd1);
2760 glVertex3fv(crd2);
2761 }
2762
2763 glEnd();
2764 glDisable(GL_LINE_STIPPLE);
2765 }
2766 }
2767 else
2768 {
2769 for (i32u n2 = 0;n2 < susf->hi_vector.size();n2++) // helix constraints
2770 {
2771 glEnable(GL_LINE_STIPPLE);
2772 glLineStipple(1, 0x3333);
2773
2774 glBegin(GL_LINES);
2775 glColor3f(0.20, 1.00, 0.10);
2776
2777 const int sz = susf->hi_vector[n2].mc_H_don.size();
2778 if (sz != (int) susf->hi_vector[n2].mc_H_acc.size()) assertion_failed(__FILE__, __LINE__, "mc_H_don/mc_H_acc mismatch");
2779
2780 for (i32u n3 = 0;n3 < sz;n3++)
2781 {
2782 const fGL * crd1 = susf->hi_vector[n2].mc_H_don[n3]->GetCRD(n1);
2783 const fGL * crd2 = susf->hi_vector[n2].mc_H_acc[n3]->GetCRD(n1);
2784
2785 glVertex3fv(crd1);
2786 glVertex3fv(crd2);
2787 }
2788
2789 glEnd();
2790 glDisable(GL_LINE_STIPPLE);
2791 }
2792
2793 for (i32u n2 = 0;n2 < susf->sp_vector.size();n2++) // strand constraints
2794 {
2795 glEnable(GL_LINE_STIPPLE);
2796 glLineStipple(1, 0x3333);
2797
2798 glBegin(GL_LINES);
2799 glColor3f(0.20, 1.00, 0.10);
2800
2801 const int sz = susf->sp_vector[n2].mc_S_don.size();
2802 if (sz != (int) susf->sp_vector[n2].mc_S_acc.size()) assertion_failed(__FILE__, __LINE__, "mc_S_don/mc_S_acc mismatch");
2803
2804 for (i32u n3 = 0;n3 < sz;n3++)
2805 {
2806 const fGL * crd1 = susf->sp_vector[n2].mc_S_don[n3]->GetCRD(n1);
2807 const fGL * crd2 = susf->sp_vector[n2].mc_S_acc[n3]->GetCRD(n1);
2808
2809 glVertex3fv(crd1);
2810 glVertex3fv(crd2);
2811 }
2812
2813 glEnd();
2814 glDisable(GL_LINE_STIPPLE);
2815 }
2816 }
2817
2818 if (susf->mode != setup1_sf::modeUA) // the Px-specific rendering starts here...
2819 {
2820 for (i32u n2 = 0;n2 < susf->chn_vector.size();n2++) // protein chains...
2821 {
2822 if (susf->mode == setup1_sf::modeP3)
2823 {
2824 for (i32s n3 = 0;n3 < ((i32s) susf->chn_vector[n2].res_vector.size()) - 1;n3++)
2825 {
2826 i32s ind1[3] = { n2, n3 + 0, 0 };
2827 i32s ind2[3] = { n2, n3 + 1, 0 };
2828
2829 const fGL * crd1 = susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->GetCRD(n1);
2830 const fGL * crd2 = susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->GetCRD(n1);
2831 const fGL * crd[2] = { crd1, crd2 };
2832
2833 fGL col1[4] = { 0.8, 0.8, 0.6, 1.0 }; // todo...
2834 fGL col2[4] = { 0.8, 0.8, 0.6, 1.0 }; // todo...
2835 const fGL * col[2] = { col1, col2 };
2836
2837 fGL rad[2] =
2838 {
2839 susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->vdwr,
2840 susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->vdwr
2841 };
2842
2843 if (wcl->render == RENDER_WIREFRAME)
2844 {
2845 glBegin(GL_LINES);
2846 glColor3fv(col1); glVertex3fv(crd1);
2847 glColor3fv(col2); glVertex3fv(crd2);
2848 glEnd();
2849 }
2850 else if (wcl->render != RENDER_NOTHING)
2851 {
2852 glEnable(GL_LIGHTING);
2853 DrawCylinder1(crd, col, rad);
2854 glDisable(GL_LIGHTING);
2855 }
2856 }
2857 }
2858
2859 for (i32u n3 = 0;n3 < susf->chn_vector[n2].res_vector.size();n3++)
2860 {
2861 const bool is_PRO = (susf->chn_vector[n2].res_vector[n3].GetSymbol() == 'P');
2862
2863 int n_sc_bonds = susf->chn_vector[n2].res_vector[n3].GetNumA() - 1;
2864 if (susf->mode == setup1_sf::modeP5)
2865 {
2866 n_sc_bonds -= 2;
2867 if (is_PRO) n_sc_bonds++;
2868 }
2869
2870 for (i32s n4 = 0;n4 < n_sc_bonds;n4++)
2871 {
2872 i32s ind1[3] = { n2, n3, NOT_DEFINED };
2873 i32s ind2[3] = { n2, n3, NOT_DEFINED };
2874
2875 if (susf->mode == setup1_sf::modeP3)
2876 {
2877 if (!n4) { ind1[2] = 0; ind2[2] = 1; }
2878 else { ind1[2] = 1; ind2[2] = 2; }
2879 }
2880 else
2881 {
2882 if (!n4) { ind1[2] = 1; ind2[2] = 3; }
2883 else
2884 {
2885 if (!is_PRO) { ind1[2] = 3; ind2[2] = 4; }
2886 else { ind1[2] = 3; ind2[2] = 0; }
2887 }
2888 }
2889
2890 const fGL * crd1 = susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->GetCRD(n1);
2891 const fGL * crd2 = susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->GetCRD(n1);
2892 const fGL * crd[2] = { crd1, crd2 };
2893
2894 fGL col1[4] = { 0.6, 0.8, 0.8, 1.0 }; // todo...
2895 fGL col2[4] = { 0.6, 0.8, 0.8, 1.0 }; // todo...
2896 const fGL * col[2] = { col1, col2 };
2897
2898 fGL rad[2] =
2899 {
2900 susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->vdwr,
2901 susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->vdwr
2902 };
2903
2904 if (wcl->render == RENDER_WIREFRAME)
2905 {
2906 glBegin(GL_LINES);
2907 glColor3fv(col1); glVertex3fv(crd1);
2908 glColor3fv(col2); glVertex3fv(crd2);
2909 glEnd();
2910 }
2911 else if (wcl->render != RENDER_NOTHING)
2912 {
2913 glEnable(GL_LIGHTING);
2914 DrawCylinder1(crd, col, rad);
2915 glDisable(GL_LIGHTING);
2916 }
2917 }
2918 }
2919 }
2920
2921 /* for (i32s n2 = 0;n2 < (i32s) susf->dsb_vector.size();n2++) // disulphide bridges.
2922 {
2923 i32s ind1[3] = { susf->dsb_vector[n2].GetChn(0), susf->dsb_vector[n2].GetRes(0), 3 };
2924 i32s ind2[3] = { susf->dsb_vector[n2].GetChn(1), susf->dsb_vector[n2].GetRes(1), 3 };
2925
2926 const fGL * crd1 = susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->GetCRD(n1);
2927 const fGL * crd2 = susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->GetCRD(n1);
2928 const fGL * crd[2] = { crd1, crd2 };
2929
2930 fGL col1[4] = { 1.0, 1.0, 0.0, 1.0 }; // todo...
2931 fGL col2[4] = { 1.0, 1.0, 0.0, 1.0 }; // todo...
2932 const fGL * col[2] = { col1, col2 };
2933
2934 fGL rad[2] =
2935 {
2936 susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->vdwr,
2937 susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->vdwr
2938 };
2939
2940 if (wcl->render == RENDER_WIREFRAME)
2941 {
2942 glBegin(GL_LINES);
2943 glColor3fv(col1); glVertex3fv(crd1);
2944 glColor3fv(col2); glVertex3fv(crd2);
2945 glEnd();
2946 }
2947 else
2948 {
2949 glEnable(GL_LIGHTING);
2950 DrawCylinder1(crd, col, rad);
2951 glDisable(GL_LIGHTING);
2952 }
2953 } */
2954 }
2955 }
2956
2957 /*//////////////////////////////////////////////////////////////////////////////////////////////////
2958 ////////////////////////////////////////////////////////////////////////////////////////////////////
2959 glEnable(GL_LIGHTING); glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, true); glBegin(GL_QUADS);
2960 // do not take the direction from any array, but calculate it using N/C/O???
2961 for (iter_bl it1 = bond_list.begin();it1 != bond_list.end();it1++) // sf peptide dipoles...
2962 {
2963 if ((* it1).sf_pbdd < -1000.0) continue;
2964
2965 //char symbol2 = chn_vector[n2].res_vector[n3 + 1].symbol;
2966 //if (symbol2 == 'P') continue; // skip all X-pro cases !!!
2967
2968 atom * prev = NULL;
2969 atom * curr = (* it1).atmr[0];
2970 atom * next = (* it1).atmr[1];
2971
2972 // WARNING!!! this is pretty slow!!! need to find the previous c-alpha.
2973 iter_cl it2;
2974 for (it2 = curr->cr_list.begin();it2 != curr->cr_list.end();it2++)
2975 {
2976 if ((* it2).atmr == next) continue;
2977
2978 if ((* it2).atmr->el.GetAtomicNumber() > 0) continue;
2979 if ((* it2).atmr->sf_atmtp & 0xFF) continue;
2980
2981 prev = (* it2).atmr;
2982 break;
2983 }
2984
2985 if (!prev) continue;
2986
2987 v3d<fGL> v1(prev->GetCRD(n1), curr->GetCRD(n1));
2988 v3d<fGL> v2(curr->GetCRD(n1), next->GetCRD(n1));
2989
2990 v3d<fGL> v3 = v1.vpr(v2); v3 = v3 * (0.075 / v3.len());
2991 v3d<fGL> v4 = v3.vpr(v2); v4 = v4 * (0.075 / v4.len());
2992
2993 fGL peptide = (* it1).sf_pbdd; // this is the same for all crd_sets!!!
2994 v3d<fGL> v5 = (v3 * sin(peptide)) + (v4 * cos(peptide));
2995
2996 fGL peptnorm = peptide - M_PI / 2.0;
2997 v3d<fGL> normal = (v3 * sin(peptnorm)) + (v4 * cos(peptnorm));
2998 normal = normal / normal.len(); glNormal3fv(normal.data);
2999
3000 v3d<fGL> pvc(curr->GetCRD(n1));
3001 v3d<fGL> pv1 = pvc + (v2 * 0.5) + v5; v3d<fGL> pv2 = pvc + (v2 * 0.90);
3002 v3d<fGL> pv3 = pvc + (v2 * 0.5) - v5; v3d<fGL> pv4 = pvc + (v2 * 0.10);
3003
3004 glColor3f(1.0, 0.0, 0.0); glVertex3fv(pv1.data);
3005 glColor3f(0.0, 1.0, 0.0); glVertex3fv(pv2.data);
3006 glColor3f(0.0, 0.0, 1.0); glVertex3fv(pv3.data);
3007 glColor3f(0.0, 1.0, 0.0); glVertex3fv(pv4.data);
3008 }
3009 glEnd(); // GL_QUADS
3010 glDisable(GL_LIGHTING); glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, false);
3011 ////////////////////////////////////////////////////////////////////////////////////////////////////
3012 //////////////////////////////////////////////////////////////////////////////////////////////////*/
3013
3014 if (accum)
3015 {
3016 RenderObjects(wcl);
3017 glAccum(GL_ACCUM, cs_vector[n1]->accum_value);
3018 }
3019 }
3020 }
3021
RenderObjects(oglview_wcl * wcl)3022 void project::RenderObjects(oglview_wcl * wcl)
3023 {
3024 base_app::GetAppB()->RenderLights(wcl->GetCam());
3025
3026 for (i32u n1 = 0;n1 < object_vector.size();n1++)
3027 {
3028 if (object_vector[n1]->transparent) continue;
3029
3030 object_vector[n1]->Render();
3031 }
3032 }
3033
BeginClientTransformation(ogl_transformer * p1)3034 void project::BeginClientTransformation(ogl_transformer * p1)
3035 {
3036 i32s sum = 0;
3037 p1->GetLD()->crd[0] = 0.0;
3038 p1->GetLD()->crd[1] = 0.0;
3039 p1->GetLD()->crd[2] = 0.0;
3040
3041 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
3042 {
3043 if (!((* it1).flags & ATOMFLAG_USER_SELECTED)) continue;
3044 for (i32u n1 = 0;n1 < cs_vector.size();n1++)
3045 {
3046 sum++;
3047 const fGL * cdata = (* it1).GetCRD(n1);
3048 p1->GetLD()->crd[0] += cdata[0];
3049 p1->GetLD()->crd[1] += cdata[1];
3050 p1->GetLD()->crd[2] += cdata[2];
3051 }
3052 }
3053
3054 if (!sum) return;
3055
3056 p1->GetLD()->crd[0] /= (fGL) sum;
3057 p1->GetLD()->crd[1] /= (fGL) sum;
3058 p1->GetLD()->crd[2] /= (fGL) sum;
3059
3060 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
3061 {
3062 if (!((* it1).flags & ATOMFLAG_USER_SELECTED)) continue;
3063 for (i32u n1 = 0;n1 < cs_vector.size();n1++)
3064 {
3065 const fGL * cdata = (* it1).GetCRD(n1);
3066
3067 fGL x = cdata[0] - p1->GetSafeLD()->crd[0];
3068 fGL y = cdata[1] - p1->GetSafeLD()->crd[1];
3069 fGL z = cdata[2] - p1->GetSafeLD()->crd[2];
3070
3071 (* it1).SetCRD(n1, x, y, z);
3072 }
3073 }
3074 }
3075
EndClientTransformation(ogl_transformer * p1)3076 void project::EndClientTransformation(ogl_transformer * p1)
3077 {
3078 fGL matrix[16]; p1->GetMatrix(matrix);
3079
3080 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
3081 {
3082 if (!((* it1).flags & ATOMFLAG_USER_SELECTED)) continue;
3083
3084 for (i32u n1 = 0;n1 < cs_vector.size();n1++)
3085 {
3086 oglv3d<fGL> posv = oglv3d<fGL>((* it1).GetCRD(n1));
3087 TransformVector(posv, matrix);
3088
3089 (* it1).SetCRD(n1, posv[0], posv[1], posv[2]);
3090 }
3091 }
3092
3093 UpdateAllGraphicsViews(); // re-draw the bonds across selection boundary!!!
3094 }
3095
DrawEvent(oglview_wcl * oglwcl,vector<iGLu> & names)3096 void project::DrawEvent(oglview_wcl * oglwcl, vector<iGLu> & names)
3097 {
3098 if (mouseinfo::button == mouseinfo::bRight) return; // the right button is for popup menus...
3099 if (project::background_job_running) return; // protect the model-data during background jobs...
3100
3101 i32s mouse[2] =
3102 {
3103 mouseinfo::latest_x,
3104 mouseinfo::latest_y
3105 };
3106
3107 if (mouseinfo::state == mouseinfo::sDown)
3108 {
3109 if (names.size() > 1 && names[0] == GLNAME_MD_TYPE1)
3110 {
3111 draw_data[0] = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
3112 }
3113 else
3114 {
3115 fGL tmp1[3]; oglwcl->GetCRD(mouse, tmp1);
3116 atom newatom(element::current_element, tmp1, cs_vector.size());
3117
3118 AddAtom_lg(newatom);
3119
3120 draw_data[0] = & atom_list.back();
3121 }
3122 }
3123 else
3124 {
3125 if (names.size() > 1 && names[0] == GLNAME_MD_TYPE1)
3126 {
3127 draw_data[1] = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
3128 }
3129 else
3130 {
3131 fGL tmp1[3]; oglwcl->GetCRD(mouse, tmp1);
3132 atom newatom(element::current_element, tmp1, cs_vector.size());
3133
3134 AddAtom_lg(newatom);
3135
3136 draw_data[1] = & atom_list.back();
3137 }
3138
3139 // if different: update bondtype or add a new bond.
3140 // if not different: change atom to different element.
3141
3142 if (draw_data[0] != draw_data[1])
3143 {
3144 bond newbond(draw_data[0], draw_data[1], bondtype::current_bondtype);
3145 iter_bl it1 = find(bond_list.begin(), bond_list.end(), newbond);
3146 if (it1 != bond_list.end())
3147 {
3148 SystemWasModified();
3149
3150 (* it1).bt = bondtype::current_bondtype;
3151
3152 custom_app::GetAppC()->BondUpdateItem(& (* it1));
3153 }
3154 else AddBond(newbond);
3155 }
3156 else
3157 {
3158 SystemWasModified();
3159
3160 draw_data[0]->el = element::current_element;
3161 draw_data[0]->mass = element::current_element.GetAtomicMass(); // also need to update these...
3162 draw_data[0]->vdwr = element::current_element.GetVDWRadius(); // also need to update these...
3163
3164 custom_app::GetAppC()->AtomUpdateItem(draw_data[0]);
3165 }
3166
3167 UpdateAllGraphicsViews();
3168 }
3169 }
3170
EraseEvent(oglview_wcl * oglwcl,vector<iGLu> & names)3171 void project::EraseEvent(oglview_wcl * oglwcl, vector<iGLu> & names)
3172 {
3173 if (mouseinfo::button == mouseinfo::bRight) return; // the right button is for popup menus...
3174 if (project::background_job_running) return; // protect the model-data during background jobs...
3175
3176 if (mouseinfo::state == mouseinfo::sDown)
3177 {
3178 if (names.size() > 1 && names[0] == GLNAME_MD_TYPE1)
3179 {
3180 draw_data[0] = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
3181 }
3182 else
3183 {
3184 draw_data[0] = NULL;
3185 }
3186 }
3187 else
3188 {
3189 if (names.size() > 1 && names[0] == GLNAME_MD_TYPE1)
3190 {
3191 draw_data[1] = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
3192 }
3193 else
3194 {
3195 draw_data[1] = NULL;
3196 }
3197
3198 if (!draw_data[0] || !draw_data[1]) return;
3199
3200 // if different: try to find and remove a bond ; may or may not succeed.
3201 // if not different: try to find and remove an atom ; SHOULD ALWAYS SUCCEED!
3202
3203 if (draw_data[0] != draw_data[1])
3204 {
3205 bond tmpbond(draw_data[0], draw_data[1], bondtype::current_bondtype);
3206 iter_bl it1 = find(bond_list.begin(), bond_list.end(), tmpbond);
3207
3208 if (it1 != bond_list.end()) RemoveBond(it1);
3209 else return;
3210 }
3211 else
3212 {
3213 iter_al it1 = find(atom_list.begin(), atom_list.end(), (* draw_data[0]));
3214
3215 if (it1 != atom_list.end())
3216 {
3217 RemoveAtom(it1);
3218
3219 // removing an atom will cause changes in atom indexing -> must update all atoms and bonds in pv!!!
3220 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3221
3222 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
3223 {
3224 custom_app::GetAppC()->AtomUpdateItem(& (* it1));
3225 }
3226
3227 for (iter_bl it1 = bond_list.begin();it1 != bond_list.end();it1++)
3228 {
3229 custom_app::GetAppC()->BondUpdateItem(& (* it1));
3230 }
3231 }
3232 else
3233 {
3234 assertion_failed(__FILE__, __LINE__, "atom not found.");
3235 }
3236 }
3237
3238 UpdateAllGraphicsViews();
3239 }
3240 }
3241
SelectEvent(oglview_wcl *,vector<iGLu> & names)3242 void project::SelectEvent(oglview_wcl *, vector<iGLu> & names)
3243 {
3244 if (project::background_job_running) return; // protect the model-data during background jobs...
3245
3246 if (names[0] == GLNAME_MD_TYPE1)
3247 {
3248 atom * ref = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
3249 bool selected = (ref->flags & ATOMFLAG_USER_SELECTED);
3250
3251 if (custom_app::GetCurrentSelectMode() == custom_app::smResidue || custom_app::GetCurrentSelectMode() == custom_app::smChain)
3252 {
3253 if (ref_civ == NULL)
3254 {
3255 ostringstream str;
3256 str << _("Group information about chains/residues is needed for this operation.") << endl;
3257 str << _("Is it OK to update group information?") << ends;
3258
3259 bool update = Question(str.str().c_str());
3260 if (!update) return;
3261
3262 UpdateChains();
3263 }
3264
3265 bool no_info = false;
3266 if (ref->id[1] == NOT_DEFINED) no_info = true;
3267 if (custom_app::GetCurrentSelectMode() == custom_app::smResidue && ref->id[2] == NOT_DEFINED) no_info = true;
3268
3269 if (no_info)
3270 {
3271 Message(_("Sorry, no chain/residue information available for this atom."));
3272 return;
3273 }
3274 }
3275
3276 if (custom_app::GetCurrentSelectMode() == custom_app::smMolecule)
3277 {
3278 if (!IsGroupsClean()) UpdateGroups();
3279 }
3280
3281 iter_al it1;
3282 iter_al range1[2];
3283 iter_al range2[2];
3284
3285 if (selected) cout << "de";
3286 switch (custom_app::GetCurrentSelectMode())
3287 {
3288 case custom_app::smAtom:
3289 ref->flags ^= ATOMFLAG_USER_SELECTED;
3290 cout << _("selected atom ") << ref->index << _(" (atomtype = ") << hex << ref->atmtp << dec << ")." << endl;
3291 break;
3292
3293 case custom_app::smResidue:
3294 GetRange(1, ref->id[1], range1); // get the chain!!!
3295 GetRange(2, range1, ref->id[2], range2); // get the residue!!!
3296 for (it1 = range2[0];it1 != range2[1];it1++)
3297 {
3298 if (selected) (* it1).flags &= (~ATOMFLAG_USER_SELECTED);
3299 else (* it1).flags |= (ATOMFLAG_USER_SELECTED);
3300 }
3301
3302 cout << _("selected residue ") << ref->id[2] << _(" from chain ") << ref->id[1] << "." << endl;
3303 break;
3304
3305 case custom_app::smChain:
3306 GetRange(1, ref->id[1], range1); // get the chain!!!
3307 for (it1 = range1[0];it1 != range1[1];it1++)
3308 {
3309 if (selected) (* it1).flags &= (~ATOMFLAG_USER_SELECTED);
3310 else (* it1).flags |= (ATOMFLAG_USER_SELECTED);
3311 }
3312
3313 cout << _("selected chain ") << ref->id[1] << "." << endl;
3314 break;
3315
3316 case custom_app::smMolecule:
3317 if (IsGroupsSorted()) // if atom_list is sorted, a quicker method based on model::GetRange() is used.
3318 {
3319 GetRange(0, ref->id[0], range1); // get the molecule!!!
3320 for (it1 = range1[0];it1 != range1[1];it1++)
3321 {
3322 if (selected) (* it1).flags &= (~ATOMFLAG_USER_SELECTED);
3323 else (* it1).flags |= (ATOMFLAG_USER_SELECTED);
3324 }
3325 }
3326 else
3327 {
3328 for (it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
3329 {
3330 if ((* it1).id[0] != ref->id[0]) continue;
3331
3332 if (selected) (* it1).flags &= (~ATOMFLAG_USER_SELECTED);
3333 else (* it1).flags |= (ATOMFLAG_USER_SELECTED);
3334 }
3335 }
3336
3337 cout << _("selected molecule ") << ref->id[0] << "." << endl;
3338 break;
3339 }
3340
3341 UpdateAllGraphicsViews();
3342 }
3343 }
3344
MeasureEvent(oglview_wcl *,vector<iGLu> & names)3345 void project::MeasureEvent(oglview_wcl *, vector<iGLu> & names)
3346 {
3347 if (project::background_job_running) return; // protect the model-data during background jobs...
3348
3349 ostringstream str1;
3350
3351 // PLEASE NOTE!!! we use always the 1st coordinate set here...
3352 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3353
3354 // we can be sure that "ref" is always up-to-date but the stored pointers
3355 // mt_a1/2/3 can be invalid ; so check them before use. reset if problems.
3356
3357 if (names[0] == GLNAME_MD_TYPE1)
3358 {
3359 atom * ref = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
3360 ref->flags ^= ATOMFLAG_MEASURE_TOOL_SEL;
3361 UpdateAllGraphicsViews();
3362
3363 if (mt_a1 == NULL)
3364 {
3365 mt_a1 = ref; // this must be OK.
3366 str1 << _("charge: ") << ref->charge << endl << ends;
3367 PrintToLog(str1.str().c_str());
3368 }
3369 else if (mt_a1 != NULL && mt_a2 == NULL)
3370 {
3371 if (mt_a1 == ref) { mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a1 = NULL; return; }
3372
3373 mt_a2 = ref; // this must be OK.
3374
3375 iter_al itX;
3376 for (itX = atom_list.begin();itX != atom_list.end();itX++) if (& (* itX) == mt_a1) break; if (itX == atom_list.end()) goto reset_all;
3377
3378 const fGL * p1 = mt_a1->GetCRD(0);
3379 const fGL * p2 = mt_a2->GetCRD(0);
3380
3381 v3d<fGL> v1(p1, p2);
3382 fGL len = v1.len();
3383
3384 str1 << _("distance: ") << len << " nm" << endl << ends;
3385 PrintToLog(str1.str().c_str());
3386 }
3387 else if (mt_a1 != NULL && mt_a2 != NULL && mt_a3 == NULL)
3388 {
3389 if (mt_a1 == ref) { mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a1 = mt_a2; mt_a2 = NULL; return; }
3390 else if (mt_a2 == ref) { mt_a2->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a2 = NULL; return; }
3391
3392 mt_a3 = ref; // this must be OK.
3393
3394 iter_al itX;
3395 for (itX = atom_list.begin();itX != atom_list.end();itX++) if (& (* itX) == mt_a1) break; if (itX == atom_list.end()) goto reset_all;
3396 for (itX = atom_list.begin();itX != atom_list.end();itX++) if (& (* itX) == mt_a2) break; if (itX == atom_list.end()) goto reset_all;
3397
3398 const fGL * p1 = mt_a1->GetCRD(0);
3399 const fGL * p2 = mt_a2->GetCRD(0);
3400 const fGL * p3 = mt_a3->GetCRD(0);
3401
3402 v3d<fGL> v1(p2, p1);
3403 v3d<fGL> v2(p2, p3);
3404 fGL ang = v1.ang(v2) * 180.0 / M_PI;
3405
3406 str1 << _("angle: ") << ang << _(" deg") << endl << ends;
3407 PrintToLog(str1.str().c_str());
3408 }
3409 else
3410 {
3411 if (mt_a1 == ref) { mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a1 = mt_a2; mt_a2 = mt_a3; mt_a3 = NULL; return; }
3412 else if (mt_a2 == ref) { mt_a2->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a2 = mt_a3; mt_a3 = NULL; return; }
3413 else if (mt_a3 == ref) { mt_a3->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a3 = NULL; return; }
3414
3415 const fGL * p1 = mt_a1->GetCRD(0);
3416 const fGL * p2 = mt_a2->GetCRD(0);
3417 const fGL * p3 = mt_a3->GetCRD(0);
3418 const fGL * p4 = ref->GetCRD(0);
3419
3420 v3d<fGL> v1(p2, p1);
3421 v3d<fGL> v2(p2, p3);
3422 v3d<fGL> v3(p3, p4);
3423 fGL tor = v1.tor(v2, v3) * 180.0 / M_PI;
3424
3425 str1 << _("torsion: ") << tor << _(" deg ") << endl << ends;
3426 PrintToLog(str1.str().c_str());
3427
3428 mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
3429 mt_a2->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
3430 mt_a3->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
3431 ref->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
3432
3433 goto reset_all;
3434 }
3435
3436 return;
3437 }
3438
3439 reset_all:
3440
3441 mt_a1 = mt_a2 = mt_a3 = NULL;
3442 UpdateAllGraphicsViews();
3443 }
3444
DoFormula(void)3445 void project::DoFormula(void)
3446 {
3447 int i;
3448 double molweight = 0.0;
3449
3450 ostringstream str;
3451
3452 int count[ELEMENT_SYMBOLS];
3453
3454 // These are the atomic numbers of the elements in alphabetical order.
3455 const int alphabetical[ELEMENT_SYMBOLS] = {
3456 89, 47, 13, 95, 18, 33, 85, 79, 5, 56, 4, 107, 83, 97, 35, 6, 20, 48,
3457 58, 98, 17, 96, 27, 24, 55, 29, 105, 66, 68, 99, 63, 9, 26, 100, 87, 31,
3458 64, 32, 1, 2, 72, 80, 67, 108, 53, 49, 77, 19, 36, 57, 3, 103, 71, 101,
3459 12, 25, 42, 109, 7, 11, 41, 60, 10, 28, 102, 93, 8, 76, 15, 91, 82, 46,
3460 61, 84, 59, 78, 94, 88, 37, 75, 104, 45, 86, 44, 16, 51, 21, 34, 106, 14,
3461 62, 50, 38, 73, 65, 43, 52, 90, 22, 81, 69, 92, 110, 23, 74, 54, 39, 70,
3462 30, 40
3463 };
3464
3465 int index;
3466
3467 for (i = 0;i < ELEMENT_SYMBOLS;i++)
3468 {
3469 count[i] = 0;
3470 }
3471
3472 iter_al it2 = atom_list.begin();
3473 while (it2 != atom_list.end())
3474 {
3475 iter_al it3 = it2++;
3476 count[(* it3).el.GetAtomicNumber() - 1]++;
3477 molweight += (* it3).mass;
3478 }
3479
3480 for (i = 0;i < ELEMENT_SYMBOLS;i++)
3481 {
3482 index = alphabetical[i] - 1;
3483 if (count[index] > 1)
3484 {
3485 str << (element(index + 1)).GetSymbol() << count[index] << " ";
3486 }
3487 else if (count[index] == 1)
3488 {
3489 str << (element(index + 1)).GetSymbol();
3490 }
3491 }
3492
3493 str << endl;
3494 str << _("MW: ") << molweight << ends;
3495
3496 Message(str.str().c_str());
3497 }
3498
DoEnergyPlot1D(i32s inda,i32s indb,i32s indc,i32s indd,i32s div1,fGL start1,fGL end1,i32s optsteps)3499 void project::DoEnergyPlot1D(i32s inda, i32s indb, i32s indc, i32s indd, i32s div1, fGL start1, fGL end1, i32s optsteps)
3500 {
3501 // 2003-11-17 : for IC modification and structure
3502 // refinement, make a temporary molecular mechanics model.
3503
3504 // 2007-01-15 : also make SF setups work (for debugging purposes).
3505
3506 setup * tmpsu = GetCurrentSetup();
3507 setup1_mm * tmpsuMM = dynamic_cast<setup1_mm *>(tmpsu);
3508 setup1_sf * tmpsuSF = dynamic_cast<setup1_sf *>(tmpsu);
3509
3510 // if current setup is not a QM one, get the eng class...
3511
3512 i32s curr_eng_index = 0;
3513 if (tmpsuMM != NULL) curr_eng_index = GetCurrentSetup()->GetCurrEngIndex();
3514 if (tmpsuSF != NULL) curr_eng_index = GetCurrentSetup()->GetCurrEngIndex();
3515
3516 model * tmpmdl = new model(); // the default setup here will be molecular mechanics!
3517
3518 vector<atom *> av; vector<atom *> av_tmp;
3519
3520 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
3521 {
3522 atom newatm((* it1).el, (* it1).GetCRD(0), tmpmdl->GetCRDSetCount());
3523 tmpmdl->AddAtom_lg(newatm);
3524
3525 av.push_back(& (* it1));
3526 av_tmp.push_back(& tmpmdl->GetLastAtom());
3527 }
3528
3529 for (iter_bl it1 = GetBondsBegin();it1 != GetBondsEnd();it1++)
3530 {
3531 i32u ind1 = 0;
3532 while (ind1 < av.size() && av[ind1] != (* it1).atmr[0]) ind1++;
3533 if (ind1 == av.size()) assertion_failed(__FILE__, __LINE__, "atom #1 not found.");
3534
3535 i32u ind2 = 0;
3536 while (ind2 < av.size() && av[ind2] != (* it1).atmr[1]) ind2++;
3537 if (ind2 == av.size()) assertion_failed(__FILE__, __LINE__, "atom #2 not found.");
3538
3539 bond newbnd(av_tmp[ind1], av_tmp[ind2], (* it1).bt);
3540 tmpmdl->AddBond(newbnd);
3541 }
3542
3543 if (tmpsuSF != NULL)
3544 {
3545 // this is for SF only ; CHECK THIS ; MIGHT BE OBSOLETE...
3546 tmpmdl->ReplaceCurrentSetup(new setup1_sf(tmpmdl, tmpsuSF->mode, false));
3547 }
3548
3549 engine * tmpeng = tmpmdl->GetCurrentSetup()->CreateEngineByIndex(curr_eng_index);
3550
3551 // the temporary model is now ok, continue...
3552
3553 engine * eng = GetCurrentSetup()->GetCurrentEngine();
3554 if (!eng) eng = GetCurrentSetup()->CreateEngineByIndex(GetCurrentSetup()->GetCurrEngIndex());
3555
3556 i32s molnum = 0; i32s in_crdset = 0;
3557
3558 i32s atmi1[4] = { inda, indb, indc, indd };
3559 atom * atmr1[4]; f64 range1[2];
3560 range1[0] = M_PI * start1 / 180.0;
3561 range1[1] = M_PI * end1 / 180.0;
3562
3563 for (i32s n1 = 0;n1 < 4;n1++)
3564 {
3565 iter_al it1;
3566
3567 it1 = tmpmdl->FindAtomByIndex(atmi1[n1]);
3568 if (it1 == tmpmdl->GetAtomsEnd())
3569 {
3570 ostringstream strE;
3571 strE << _("ERROR : atom ") << (n1 + 1) << _(" not found!") << ends;
3572
3573 PrintToLog(strE.str().c_str());
3574 return;
3575 }
3576
3577 atmr1[n1] = & (* it1);
3578 }
3579
3580 // must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
3581 // must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
3582 // must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
3583 if (!tmpmdl->IsGroupsClean()) tmpmdl->UpdateGroups(); // for internal coordinates...
3584 if (!tmpmdl->IsGroupsSorted()) tmpmdl->SortGroups(); // for internal coordinates...
3585
3586 intcrd * tmpic = new intcrd((* tmpmdl), molnum, in_crdset);
3587
3588 i32s ict1 = tmpic->FindTorsion(atmr1[1], atmr1[2]);
3589 if (ict1 < 0 && tmpsuSF == NULL)
3590 {
3591 PrintToLog(_("ERROR : could not find ic.\n"));
3592 return;
3593 }
3594
3595 if (tmpsuSF != NULL)
3596 {
3597 // this is for SF only...
3598 CopyCRD(tmpmdl, tmpeng, 0);
3599 }
3600
3601 v3d<fGL> v1a(atmr1[1]->GetCRD(in_crdset), atmr1[0]->GetCRD(in_crdset));
3602 v3d<fGL> v1b(atmr1[1]->GetCRD(in_crdset), atmr1[2]->GetCRD(in_crdset));
3603 v3d<fGL> v1c(atmr1[2]->GetCRD(in_crdset), atmr1[3]->GetCRD(in_crdset));
3604 f64 oldt1 = v1a.tor(v1b, v1c);
3605
3606 bool success = tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], oldt1, 10.0, false);
3607 if (!success)
3608 {
3609 PrintToLog(_("ERROR : could not find tor-term.\n"));
3610 return;
3611 }
3612
3613 const char * s1 = _("tor(deg)"); const char * sv = _("E(kJ/mol)");
3614 p1dview_wcl * plot = AddPlot1DClient(s1, sv, true);
3615
3616 f64 tor1 = range1[0];
3617 for (i32s s1 = 0;s1 < (div1 + 1);s1++)
3618 {
3619 if (ict1 < 0)
3620 {
3621 // this is for SF only...
3622 tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], tor1, 5000.0, true);
3623 }
3624 else
3625 {
3626 tmpic->SetTorsion(ict1, tor1 - oldt1);
3627 tmpic->UpdateCartesian();
3628
3629 CopyCRD(tmpmdl, tmpeng, 0); // lock_local_structure needs coordinates!!!
3630 tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], tor1, 5000.0, true);
3631 }
3632
3633 // optimize...
3634
3635 geomopt * opt = new geomopt(tmpeng, 100, 0.025, 10.0); // optimal settings?!?!?
3636
3637 for (i32s n1 = 0;n1 < optsteps;n1++)
3638 {
3639 opt->TakeCGStep(conjugate_gradient::Newton2An);
3640 if (!(n1 % 50)) cout << n1 << " " << opt->optval << " " << opt->optstp << endl;
3641 }
3642
3643 CopyCRD(tmpeng, tmpmdl, 0);
3644 tmpmdl->CenterCRDSet(0, true);
3645 delete opt;
3646
3647 for (i32u n1 = 0;n1 < av_tmp.size();n1++)
3648 {
3649 const fGL * tmpcrd = av_tmp[n1]->GetCRD(0);
3650 av[n1]->SetCRD(0, tmpcrd[0], tmpcrd[1], tmpcrd[2]);
3651 }
3652
3653 // compute energy for final structure...
3654
3655 f64 value;
3656 CopyCRD(this, eng, 0);
3657 eng->Compute(0); value = eng->GetEnergy();
3658
3659 // ...and add it to the plot.
3660
3661 plot->AddDataWithAC(180.0 * tor1 / M_PI, value, eng);
3662
3663 ostringstream str1;
3664 str1 << _("tor = ") << (180.0 * tor1 / M_PI) << _(" deg, energy = ") << value << _(" kJ/mol.") << endl << ends;
3665 PrintToLog(str1.str().c_str());
3666
3667 tor1 += (range1[1] - range1[0]) / (f64) div1;
3668 }
3669
3670 delete tmpic;
3671 delete tmpeng;
3672 delete tmpmdl;
3673
3674 // the "eng" object is the setup->current_eng object, so there's no need to delete it...
3675
3676 plot->Finalize();
3677 plot->GetWnd()->RequestUpdate(false);
3678 }
3679
DoEnergyPlot2D(i32s inda,i32s indb,i32s indc,i32s indd,i32s div1,fGL start1,fGL end1,i32s indi,i32s indj,i32s indk,i32s indl,i32s div2,fGL start2,fGL end2,i32s optsteps)3680 void project::DoEnergyPlot2D(i32s inda, i32s indb, i32s indc, i32s indd, i32s div1, fGL start1, fGL end1, i32s indi, i32s indj, i32s indk, i32s indl, i32s div2, fGL start2, fGL end2, i32s optsteps)
3681 {
3682 // 2003-11-17 : for IC modification and structure
3683 // refinement, make a temporary molecular mechanics model.
3684
3685 // 2007-01-15 : also make SF setups work (for debugging purposes).
3686
3687 setup * tmpsu = GetCurrentSetup();
3688 setup1_mm * tmpsuMM = dynamic_cast<setup1_mm *>(tmpsu);
3689 setup1_sf * tmpsuSF = dynamic_cast<setup1_sf *>(tmpsu);
3690
3691 // if current setup is not a QM one, get the eng class...
3692
3693 i32s curr_eng_index = 0;
3694 if (tmpsuMM != NULL) curr_eng_index = GetCurrentSetup()->GetCurrEngIndex();
3695 if (tmpsuSF != NULL) curr_eng_index = GetCurrentSetup()->GetCurrEngIndex();
3696
3697 model * tmpmdl = new model(); // the default setup here will be molecular mechanics!
3698
3699 vector<atom *> av; vector<atom *> av_tmp;
3700
3701 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
3702 {
3703 atom newatm((* it1).el, (* it1).GetCRD(0), tmpmdl->GetCRDSetCount());
3704 tmpmdl->AddAtom_lg(newatm);
3705
3706 av.push_back(& (* it1));
3707 av_tmp.push_back(& tmpmdl->GetLastAtom());
3708 }
3709
3710 for (iter_bl it1 = GetBondsBegin();it1 != GetBondsEnd();it1++)
3711 {
3712 i32u ind1 = 0;
3713 while (ind1 < av.size() && av[ind1] != (* it1).atmr[0]) ind1++;
3714 if (ind1 == av.size()) assertion_failed(__FILE__, __LINE__, "atom #1 not found.");
3715
3716 i32u ind2 = 0;
3717 while (ind2 < av.size() && av[ind2] != (* it1).atmr[1]) ind2++;
3718 if (ind2 == av.size()) assertion_failed(__FILE__, __LINE__, "atom #2 not found.");
3719
3720 bond newbnd(av_tmp[ind1], av_tmp[ind2], (* it1).bt);
3721 tmpmdl->AddBond(newbnd);
3722 }
3723
3724 if (tmpsuSF != NULL)
3725 {
3726 // this is for SF only ; CHECK THIS ; MIGHT BE OBSOLETE...
3727 tmpmdl->ReplaceCurrentSetup(new setup1_sf(tmpmdl, tmpsuSF->mode, false));
3728 }
3729
3730 engine * tmpeng = tmpmdl->GetCurrentSetup()->CreateEngineByIndex(curr_eng_index);
3731
3732 // the temporary model is now ok, continue...
3733
3734 engine * eng = GetCurrentSetup()->GetCurrentEngine();
3735 if (!eng) eng = GetCurrentSetup()->CreateEngineByIndex(GetCurrentSetup()->GetCurrEngIndex());
3736
3737 i32s molnum = 0; i32s in_crdset = 0;
3738
3739 i32s atmi1[4] = { inda, indb, indc, indd };
3740 atom * atmr1[4]; f64 range1[2];
3741 range1[0] = M_PI * start1 / 180.0;
3742 range1[1] = M_PI * end1 / 180.0;
3743
3744 i32s atmi2[4] = { indi, indj, indk, indl };
3745 atom * atmr2[4]; f64 range2[2];
3746 range2[0] = M_PI * start2 / 180.0;
3747 range2[1] = M_PI * end2 / 180.0;
3748
3749 for (i32s n1 = 0;n1 < 4;n1++)
3750 {
3751 iter_al it1;
3752
3753 it1 = tmpmdl->FindAtomByIndex(atmi1[n1]);
3754 if (it1 == tmpmdl->GetAtomsEnd())
3755 {
3756 ostringstream strE;
3757 strE << _("ERROR : tor1 atom ") << (n1 + 1) << _(" not found!") << endl << ends;
3758
3759 PrintToLog(strE.str().c_str());
3760 return;
3761 }
3762
3763 atmr1[n1] = & (* it1);
3764
3765 it1 = tmpmdl->FindAtomByIndex(atmi2[n1]);
3766 if (it1 == tmpmdl->GetAtomsEnd())
3767 {
3768 ostringstream strE;
3769 strE << _("ERROR : tor2 atom ") << (n1 + 1) << _(" not found!") << endl << ends;
3770
3771 PrintToLog(strE.str().c_str());
3772 return;
3773 }
3774
3775 atmr2[n1] = & (* it1);
3776 }
3777
3778 // must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
3779 // must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
3780 // must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
3781 if (!tmpmdl->IsGroupsClean()) tmpmdl->UpdateGroups(); // for internal coordinates...
3782 if (!tmpmdl->IsGroupsSorted()) tmpmdl->SortGroups(); // for internal coordinates...
3783
3784 intcrd * tmpic = new intcrd((* tmpmdl), molnum, in_crdset);
3785
3786 i32s ict1 = tmpic->FindTorsion(atmr1[1], atmr1[2]);
3787 if (ict1 < 0 && tmpsuSF == NULL)
3788 {
3789 PrintToLog(_("ERROR : could not find ic for tor1.\n"));
3790 return;
3791 }
3792
3793 i32s ict2 = tmpic->FindTorsion(atmr2[1], atmr2[2]);
3794 if (ict2 < 0 && tmpsuSF == NULL)
3795 {
3796 PrintToLog(_("ERROR : could not find ic for tor2.\n"));
3797 return;
3798 }
3799
3800 if (tmpsuSF != NULL)
3801 {
3802 // this is for SF only...
3803 CopyCRD(tmpmdl, tmpeng, 0);
3804 }
3805
3806 v3d<fGL> v1a(atmr1[1]->GetCRD(in_crdset), atmr1[0]->GetCRD(in_crdset));
3807 v3d<fGL> v1b(atmr1[1]->GetCRD(in_crdset), atmr1[2]->GetCRD(in_crdset));
3808 v3d<fGL> v1c(atmr1[2]->GetCRD(in_crdset), atmr1[3]->GetCRD(in_crdset));
3809 f64 oldt1 = v1a.tor(v1b, v1c);
3810
3811 v3d<fGL> v2a(atmr2[1]->GetCRD(in_crdset), atmr2[0]->GetCRD(in_crdset));
3812 v3d<fGL> v2b(atmr2[1]->GetCRD(in_crdset), atmr2[2]->GetCRD(in_crdset));
3813 v3d<fGL> v2c(atmr2[2]->GetCRD(in_crdset), atmr2[3]->GetCRD(in_crdset));
3814 f64 oldt2 = v2a.tor(v2b, v2c);
3815
3816 bool success1 = tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], oldt1, 10.0, false);
3817 if (!success1)
3818 {
3819 PrintToLog(_("ERROR : could not find tor-term for tor1.\n"));
3820 return;
3821 }
3822
3823 bool success2 = tmpeng->SetTorsionConstraint(atmr2[0], atmr2[1], atmr2[2], atmr2[3], oldt2, 10.0, false);
3824 if (!success2)
3825 {
3826 PrintToLog(_("ERROR : could not find tor-term for tor2.\n"));
3827 return;
3828 }
3829
3830 const char * s1 = _("tor1(deg)"); const char * s2 = _("tor2(deg)"); const char * sv = _("E(kJ/mol)");
3831 p2dview_wcl * plot = AddPlot2DClient(s1, s2, sv, true);
3832
3833 f64 tor1 = range1[0];
3834 for (i32s s1 = 0;s1 < (div1 + 1);s1++)
3835 {
3836 f64 tor2 = range2[0];
3837 for (i32s s2 = 0;s2 < (div2 + 1);s2++)
3838 {
3839 if (ict1 < 0 || ict2 < 0)
3840 {
3841 // this is for SF only...
3842 tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], tor1, 5000.0, true);
3843 tmpeng->SetTorsionConstraint(atmr2[0], atmr2[1], atmr2[2], atmr2[3], tor2, 5000.0, true);
3844 }
3845 else
3846 {
3847 tmpic->SetTorsion(ict1, tor1 - oldt1);
3848 tmpic->SetTorsion(ict2, tor2 - oldt2);
3849 tmpic->UpdateCartesian();
3850
3851 CopyCRD(tmpmdl, tmpeng, 0); // lock_local_structure needs coordinates!!!
3852 tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], tor1, 5000.0, true);
3853 tmpeng->SetTorsionConstraint(atmr2[0], atmr2[1], atmr2[2], atmr2[3], tor2, 5000.0, true);
3854 }
3855
3856 // optimize...
3857
3858 geomopt * opt = new geomopt(tmpeng, 100, 0.025, 10.0); // optimal settings?!?!?
3859
3860 for (i32s n1 = 0;n1 < optsteps;n1++)
3861 {
3862 opt->TakeCGStep(conjugate_gradient::Newton2An);
3863 if (!(n1 % 50)) cout << n1 << " " << opt->optval << " " << opt->optstp << endl;
3864 }
3865
3866 CopyCRD(tmpeng, tmpmdl, 0);
3867 tmpmdl->CenterCRDSet(0, true);
3868 delete opt;
3869
3870 for (i32u n1 = 0;n1 < av_tmp.size();n1++)
3871 {
3872 const fGL * tmpcrd = av_tmp[n1]->GetCRD(0);
3873 av[n1]->SetCRD(0, tmpcrd[0], tmpcrd[1], tmpcrd[2]);
3874 }
3875
3876 // compute energy for final structure...
3877
3878 f64 value;
3879 CopyCRD(this, eng, 0);
3880 eng->Compute(0); value = eng->GetEnergy();
3881
3882 // ...and add it to the plot.
3883
3884 plot->AddDataWithAC(180.0 * tor1 / M_PI, 180.0 * tor2 / M_PI, value, eng);
3885
3886 ostringstream str1;
3887 str1 << _("tor1 = ") << (180.0 * tor1 / M_PI) << _(" deg, tor2 = ") << (180.0 * tor2 / M_PI) << _(" deg, energy = ") << value << _(" kJ/mol.") << endl << ends;
3888 PrintToLog(str1.str().c_str());
3889
3890 tor2 += (range2[1] - range2[0]) / (f64) div2;
3891 }
3892
3893 tor1 += (range1[1] - range1[0]) / (f64) div1;
3894 }
3895
3896 delete tmpic;
3897 delete tmpeng;
3898 delete tmpmdl;
3899
3900 // the "eng" object is the setup->current_eng object, so there's no need to delete it...
3901
3902 plot->Finalize();
3903 plot->GetWnd()->RequestUpdate(false);
3904 }
3905
DoTransitionStateSearch(f64 deltae,f64 initfc)3906 void project::DoTransitionStateSearch(f64 deltae, f64 initfc)
3907 {
3908 transition_state_search * tss = new transition_state_search(this, deltae, initfc);
3909 if (tss->InitFailed()) { delete tss; tss = NULL; return; }
3910
3911 ostringstream txts1;
3912 txts1 << _("r-energy = ") << tss->GetE(0) << " " << _("p-energy = ") << tss->GetE(1) << " ";
3913 txts1 << (tss->GetE(0) < tss->GetE(1) ? "r" : "p") << _(" is lower ") << fabs(tss->GetE(0) - tss->GetE(1));
3914 txts1 << endl << ends;
3915
3916 PrintToLog(txts1.str().c_str());
3917 cout << txts1.str().c_str();
3918 // char stop1; cin >> stop1;
3919
3920 f64 erl = tss->GetE(0); f64 epl = tss->GetE(1);
3921
3922 const char * s1 = "rc"; const char * sv = _("E(kJ/mol)");
3923 rcpview_wcl * plot = AddReactionCoordinatePlotClient(s1, sv, true);
3924
3925 for (i32u n1 = 0;n1 < tss->patoms.size();n1++) plot->AddPAtom(tss->patoms[n1]);
3926 for (i32u n1 = 0;n1 < tss->rbonds.size();n1++) plot->AddRBond(tss->rbonds[n1]);
3927 for (i32u n1 = 0;n1 < tss->pbonds.size();n1++) plot->AddPBond(tss->pbonds[n1]);
3928
3929 void * udata;
3930
3931 // add the initial structures...
3932 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3933
3934 plot->AddDataWithAC(tss->GetP(0), tss->GetE(0), this, 0);
3935
3936 plot->AddDataWithAC(tss->GetP(1), tss->GetE(1), this, 1);
3937
3938 // loop...
3939 // ^^^^^^^
3940
3941 i32s prev_not_stored[2] = { false, false };
3942 while (true)
3943 {
3944 tss->Run(prev_not_stored);
3945
3946 ostringstream txts2;
3947 txts2 << _("r-energy = ") << tss->GetE(0) << " " << _("p-energy = ") << tss->GetE(1) << " ";
3948 txts2 << (tss->GetE(0) < tss->GetE(1) ? "r" : "p") << _(" is lower ") << fabs(tss->GetE(0) - tss->GetE(1)) << " ";
3949
3950 if (tss->GetR(0) && tss->GetR(1))
3951 {
3952 txts2 << _("READY!") << endl << ends;
3953 PrintToLog(txts2.str().c_str());
3954 break;
3955 }
3956
3957 bool update[2] = { !tss->GetR(0), !tss->GetR(1) };
3958 if (tss->GetE(1) < erl) update[0] = false;
3959 if (tss->GetE(0) < epl) update[1] = false;
3960
3961 if (!update[0] && !update[1]) // this is a deadlock situation, fix it...
3962 {
3963 ////////////////////////////////////////////////////////////////
3964 //cout << (i32s) update[0] << (i32s) update[1] << " ";
3965 //cout << (i32s) tss->GetR(0) << (i32s) tss->GetR(1) << " ";
3966 //cout << "DEADLOCK!!!" << endl; int xx;cin>>xx;
3967 ////////////////////////////////////////////////////////////////
3968 if (!tss->GetR(0) && tss->GetR(1)) update[0] = true;
3969 if (tss->GetR(0) && !tss->GetR(1)) update[1] = true;
3970 if (!update[0] && !update[1])
3971 {
3972 f64 delta1 = erl - tss->GetE(1);
3973 f64 delta2 = epl - tss->GetE(0);
3974 i32s uuu = (delta1 > delta2 ? 0 : 1); // update the bigger one...
3975 update[uuu] = true;
3976 }
3977 }
3978
3979 txts2 << (i32s) update[0] << (i32s) update[1] << " ";
3980 txts2 << (i32s) tss->GetR(0) << (i32s) tss->GetR(1);
3981 txts2 << endl << ends;
3982
3983 PrintToLog(txts2.str().c_str());
3984 cout << txts2.str().c_str();
3985 // char stop2; cin >> stop2;
3986
3987 tss->UpdateTargets(update);
3988
3989 if (update[0])
3990 {
3991 plot->AddDataWithAC(tss->GetP(0), tss->GetE(0), this, 0);
3992 erl = tss->GetE(0);
3993 }
3994
3995 if (update[1])
3996 {
3997 plot->AddDataWithAC(tss->GetP(1), tss->GetE(1), this, 1);
3998 epl = tss->GetE(1);
3999 }
4000
4001 prev_not_stored[0] = !update[0];
4002 prev_not_stored[1] = !update[1];
4003
4004 UpdateAllGraphicsViews(true); // debug...
4005 }
4006
4007 delete tss; tss = NULL;
4008
4009 // create an approximate TS as an average of the two structures.
4010 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
4011
4012 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
4013 {
4014 const fGL * crdr = (* it1).GetCRD(0);
4015 const fGL * crdp = (* it1).GetCRD(1);
4016
4017 fGL x = (crdr[0] + crdp[0]) / 2.0;
4018 fGL y = (crdr[1] + crdp[1]) / 2.0;
4019 fGL z = (crdr[2] + crdp[2]) / 2.0;
4020
4021 (* it1).SetCRD(0, x, y, z);
4022 }
4023
4024 PopCRDSets(1); // remove the 2nd crd-set that is no longer needed.
4025
4026 // refine the approximate TS using stationary state search...
4027 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
4028
4029 DoStationaryStateSearch(100);
4030 f64 ts_e = GetCurrentSetup()->GetCurrentEngine()->GetEnergy();
4031
4032 // add the final estimate of TS, and finish the plot.
4033 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
4034
4035 plot->AddDataWithAC(0, ts_e, this, 0);
4036
4037 plot->Finalize();
4038 plot->GetWnd()->RequestUpdate(false);
4039 }
4040
DoStationaryStateSearch(i32s steps)4041 void project::DoStationaryStateSearch(i32s steps)
4042 {
4043 engine * eng = GetCurrentSetup()->GetCurrentEngine();
4044 if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
4045 eng = GetCurrentSetup()->GetCurrentEngine();
4046 if (eng == NULL) return;
4047
4048 ostringstream str1;
4049 str1 << _("Starting Stationary State Search ");
4050 str1 << _("(setup = ") << GetCurrentSetup()->GetClassName_lg();
4051 str1 << _(", engine = ") << GetCurrentSetup()->GetEngineName(GetCurrentSetup()->GetCurrEngIndex());
4052 str1 << ")." << endl << ends;
4053
4054 PrintToLog(str1.str().c_str());
4055
4056 CopyCRD(this, eng, 0);
4057
4058 // use a small default steplength; also setting maximum steplength is important!!!
4059 stationary_state_search * sss = new stationary_state_search(eng, 25, 1.0e-7, 1.0e-5);
4060
4061 char buffer[1024];
4062 PrintToLog(_("Cycle Gradient Step\n"));
4063
4064 i32s n1 = 0; // n1 counts the number of steps...
4065 while (true)
4066 {
4067 sss->TakeCGStep(conjugate_gradient::Simple);
4068
4069 sprintf(buffer, "%4d %10.4e %10.4e \n", n1, sss->optval, sss->optstp);
4070
4071 PrintToLog(buffer);
4072
4073 bool terminate = false;
4074 if (n1 >= steps)
4075 {
4076 terminate = true;
4077 PrintToLog(_("the nsteps termination test was passed.\n"));
4078 }
4079
4080 if (!(n1 % 10) || terminate)
4081 {
4082 CopyCRD(eng, this, 0);
4083 CenterCRDSet(0, true);
4084
4085 UpdateAllGraphicsViews(true);
4086 }
4087
4088 if (terminate) break; // exit the loop here!!!
4089
4090 n1++; // update the number of steps...
4091 }
4092
4093 delete sss; sss = NULL;
4094
4095 // we will not delete current_eng here, so that we can draw plots using it...
4096
4097 // above, CopyCRD was done eng->mdl and then CenterCRDSet() was done for mdl.
4098 // this might cause that old coordinates remain in eng object, possibly affecting plots.
4099 // here we sync the coordinates and other plotting data in the eng object.
4100
4101 CopyCRD(this, eng, 0);
4102 SetupPlotting();
4103 }
4104
4105 /*################################################################################################*/
4106
dummy_project(void)4107 dummy_project::dummy_project(void) :
4108 project()
4109 {
4110 }
4111
~dummy_project(void)4112 dummy_project::~dummy_project(void)
4113 {
4114 }
4115
4116 /*################################################################################################*/
4117
GetColor4(const void * dd,i32s cs,fGL * pp)4118 void color_mode_element::GetColor4(const void * dd, i32s cs, fGL * pp)
4119 {
4120 atom * ref = (atom *) dd;
4121 const fGL * color = ref->el.GetColor();
4122 pp[0] = color[0]; pp[1] = color[1]; pp[2] = color[2]; pp[3] = 1.0;
4123 }
4124
GetColor4(const void * dd,i32s cs,fGL * pp)4125 void color_mode_secstruct::GetColor4(const void * dd, i32s cs, fGL * pp)
4126 {
4127 atom * ref = (atom *) dd;
4128 model * mdl = ref->GetModel();
4129
4130 pp[0] = 0.0; pp[1] = 0.0; pp[2] = 1.0; pp[3] = 0; // loop
4131
4132 if (mdl == NULL || mdl->GetCI() == NULL) return;
4133 if (ref->id[1] < 0 || ref->id[2] < 0) return;
4134
4135 vector<chn_info> & ci_vector = (* mdl->GetCI());
4136 const char * tmptab = ci_vector[ref->id[1]].GetSecStrStates();
4137
4138 if (tmptab == NULL) return;
4139 char state = tmptab[ref->id[2]];
4140
4141 switch (state)
4142 {
4143 case '4':
4144 pp[0] = 1.0; pp[1] = 0.0; pp[2] = 0.0; // helix
4145 return;
4146
4147 case 'S':
4148 pp[0] = 0.0; pp[1] = 1.0; pp[2] = 0.0; // strand
4149 return;
4150 }
4151 }
4152
GetColor4(const void * dd,i32s cs,fGL * pp)4153 void color_mode_hydphob::GetColor4(const void * dd, i32s cs, fGL * pp)
4154 {
4155 atom * ref = (atom *) dd;
4156 model * mdl = ref->GetModel();
4157
4158 pp[0] = 0.0; pp[1] = 0.5; pp[2] = 0.0; pp[3] = 0; // default...
4159
4160 if (mdl == NULL || mdl->GetCI() == NULL) return;
4161 if (ref->id[1] < 0 || ref->id[2] < 0) return;
4162
4163 vector<chn_info> & ci_vector = (* mdl->GetCI());
4164 const char * tmp_seq1 = ci_vector[ref->id[1]].GetSequence1();
4165
4166 if (tmp_seq1 == NULL) return;
4167 char res1 = tmp_seq1[ref->id[2]];
4168
4169 switch (res1)
4170 {
4171 case 'A':
4172 case 'G':
4173 pp[0] = 0.0; pp[1] = 1.0; pp[2] = 0.0; // ala/gly
4174 return;
4175
4176 case 'V':
4177 case 'F':
4178 case 'I':
4179 case 'L':
4180 case 'P':
4181 case 'M':
4182 pp[0] = 1.0; pp[1] = 0.0; pp[2] = 0.0; // hydrophobic
4183 return;
4184
4185 case 'D':
4186 case 'E':
4187 case 'K':
4188 case 'R':
4189 pp[0] = 0.2; pp[1] = 0.2; pp[2] = 1.0; // charged
4190 return;
4191
4192 case 'S':
4193 case 'T':
4194 case 'Y':
4195 case 'C':
4196 case 'N':
4197 case 'Q':
4198 case 'H':
4199 case 'W':
4200 pp[0] = 0.0; pp[1] = 1.0; pp[2] = 2.0; // polar
4201 return;
4202 }
4203 }
4204
4205 /*################################################################################################*/
4206
4207 // eof
4208