FOX/ObjCryst++  2022
ZScatterer.cpp
1 /* ObjCryst++ Object-Oriented Crystallographic Library
2  (c) 2000-2002 Vincent Favre-Nicolin vincefn@users.sourceforge.net
3  2000-2001 University of Geneva (Switzerland)
4 
5  This program 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 program 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 program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 */
19 /* Atom.h
20 * header file for the Atom scatterer
21 *
22 */
23 
24 #include <cmath>
25 #include <typeinfo>
26 #include <cstdlib>
27 #include <stdio.h> //for sprintf()
28 
29 //#include "ObjCryst/ObjCryst/ScatteringPower.h"
30 //#include "ObjCryst/ObjCryst/Scatterer.h"
31 //#include "ObjCryst/ObjCryst/Atom.h"
32 #include "ObjCryst/ObjCryst/ZScatterer.h"
33 #include "ObjCryst/ObjCryst/ScatteringData.h"
34 
35 #include "ObjCryst/Quirks/VFNStreamFormat.h" //simple formatting of integers, REALs..
36 
37 #include "ObjCryst/Quirks/VFNDebug.h"
38 
39 #include <fstream>
40 #include <iomanip>
41 
42 #ifdef OBJCRYST_GL
43  #ifdef __DARWIN__
44  #include <OpenGL/glu.h>
45  #else
46  #include <GL/glu.h>
47  #endif
48 #endif
49 
50 #ifdef __WX__CRYST__
51  #include "ObjCryst/wxCryst/wxZScatterer.h"
52  #undef GetClassName // Conflict from wxMSW headers ? (cygwin)
53 #endif
54 
55 #ifdef _MSC_VER // MS VC++ predefined macros....
56 #undef min
57 #undef max
58 #endif
59 
60 namespace ObjCryst
61 {
62 //######################################################################
63 //
64 // ZAtom
65 //
66 //
67 //######################################################################
68 ZAtom::ZAtom(ZScatterer &scatt,const ScatteringPower *pow,
69  const long atomBond, const REAL bondLength,
70  const long atomAngle, const REAL bondAngle,
71  const long atomDihedral, const REAL dihedralAngle,
72  const REAL popu, const string &name):
73 mpScattPow(pow),
74 mAtomBond(atomBond),mAtomAngle(atomAngle),mAtomDihed(atomDihedral),
75 mBondLength(bondLength),mAngle(bondAngle),mDihed(dihedralAngle),
76 mOccupancy(popu),mName(name),mpScatt(&scatt)
77 {
78  VFN_DEBUG_MESSAGE("ZAtom::ZAtom():("<<mName<<")",5)
79 }
80 ZAtom::~ZAtom()
81 {
82  VFN_DEBUG_MESSAGE("ZAtom::~ZAtom():("<<mName<<")",5)
83 }
84 const string& ZAtom::GetClassName()const
85 {
86  static string className="ZAtom";
87  return className;
88 }
89 const string& ZAtom::GetName()const {return mName;}
90 
91 ZScatterer& ZAtom::GetZScatterer(){return *mpScatt;}
92 const ZScatterer& ZAtom::GetZScatterer()const{return *mpScatt;}
93 
94 void ZAtom::SetName(const string& name) {mName=name;}
95 long ZAtom::GetZBondAtom()const {return mAtomBond;}
96 long ZAtom::GetZAngleAtom()const {return mAtomAngle;}
97 long ZAtom::GetZDihedralAngleAtom()const {return mAtomDihed;}
98 const REAL& ZAtom::GetZBondLength()const {return mBondLength;}
99 const REAL& ZAtom::GetZAngle()const {return mAngle;}
100 const REAL& ZAtom::GetZDihedralAngle()const {return mDihed;}
101 const REAL& ZAtom::GetOccupancy()const {return mOccupancy;}
102 const ScatteringPower* ZAtom::GetScatteringPower()const{return mpScattPow;}
103 //:TODO: fix the following so that their clocks are clicked accordingly
104 void ZAtom::SetZBondLength(const REAL bond) {mBondLength=bond;mpScatt->GetClockScatterer().Click();}
105 void ZAtom::SetZAngle(const REAL angle) {mAngle=angle;mpScatt->GetClockScatterer().Click();}
106 void ZAtom::SetZDihedralAngle(const REAL dihed) {mDihed=dihed;mpScatt->GetClockScatterer().Click();}
107 void ZAtom::SetOccupancy(const REAL pop) {mOccupancy=pop;mpScatt->GetClockScatterer().Click();}
108 void ZAtom::SetScatteringPower(const ScatteringPower* scatt) {mpScattPow=scatt;mpScatt->GetClockScatterer().Click();}
109 #ifdef __WX__CRYST__
110 WXCrystObjBasic* ZAtom::WXCreate(wxWindow *parent)
111 {
112  VFN_DEBUG_ENTRY("ZAtom::WXCreate()",7)
113  mpWXCrystObj=new WXZAtom (parent,this);
114  VFN_DEBUG_EXIT("ZAtom::WXCreate()",7)
115  return mpWXCrystObj;
116 }
117 WXCrystObjBasic* ZAtom::WXGet()
118 {
119  return mpWXCrystObj;
120 }
121 void ZAtom::WXDelete()
122 {
123  if(0!=mpWXCrystObj)
124  {
125  VFN_DEBUG_MESSAGE("ZAtom::WXDelete()",5)
126  delete mpWXCrystObj;
127  }
128  mpWXCrystObj=0;
129 }
130 void ZAtom::WXNotifyDelete()
131 {
132  VFN_DEBUG_MESSAGE("ZAtom::WXNotifyDelete():"<<mName,10)
133  mpWXCrystObj=0;
134 }
135 #endif
136 //######################################################################
137 //
138 // ZMoveMinimizer
139 //
140 //
141 //######################################################################
142 ZMoveMinimizer::ZMoveMinimizer(ZScatterer &scatt):
143 RefinableObj(true),
144 mpZScatt(&scatt),mOptimObj(true)
145 {
146  mOptimObj.SetAlgorithmSimulAnnealing(ANNEALING_EXPONENTIAL,.1,.001,
147  ANNEALING_EXPONENTIAL,16,.25);
148  mOptimObj.AddRefinableObj(*this);
149 }
150 ZMoveMinimizer::~ZMoveMinimizer(){}
151 
152 
153 REAL ZMoveMinimizer::GetLogLikelihood() const
154 {
155  TAU_PROFILE("ZMoveMinimizer::GetLogLikelihood()","void ()",TAU_DEFAULT);
156  const REAL *pX1=mpZScatt->GetXCoord().data();
157  const REAL *pY1=mpZScatt->GetYCoord().data();
158  const REAL *pZ1=mpZScatt->GetZCoord().data();
159  const REAL *pX0=mXCoord0.data();
160  const REAL *pY0=mYCoord0.data();
161  const REAL *pZ0=mZCoord0.data();
162  const REAL *pW=mAtomWeight.data();
163  REAL dist=0;
164  //for(int i=mXCoord0.numElements()-1;i>=0;i--)
165  // dist+= abs(*pX1++ - *pX0++) + abs(*pY1++ - *pY0++) + abs(*pZ1++ - *pZ0++);
166  for(int i=mXCoord0.numElements()-1;i>=0;i--)
167  {
168  dist+= *pW++* ( (*pX1 - *pX0)*(*pX1 - *pX0)
169  +(*pY1 - *pY0)*(*pY1 - *pY0)
170  +(*pZ1 - *pZ0)*(*pZ1 - *pZ0));
171  pX1++;pY1++;pZ1++;
172  pX0++;pY0++;pZ0++;
173  }
174 
175  #if 0
176  const CrystVector_REAL *pXcoord=&(mpZScatt->GetXCoord());
177  const CrystVector_REAL *pYcoord=&(mpZScatt->GetYCoord());
178  const CrystVector_REAL *pZcoord=&(mpZScatt->GetZCoord());
179  REAL dist=0;
180  for(int i=pXcoord->numElements()-1;i>=0;i--)
181  {
182  dist+=mAtomWeight(i)*( ((*pXcoord)(i)-mXCoord0(i))*((*pXcoord)(i)-mXCoord0(i))
183  +((*pYcoord)(i)-mYCoord0(i))*((*pYcoord)(i)-mYCoord0(i))
184  +((*pZcoord)(i)-mZCoord0(i))*((*pZcoord)(i)-mZCoord0(i)));
185  }
186  #endif
187  return dist/mAtomWeight.sum();
188 }
189 void ZMoveMinimizer::RecordConformation()
190 {
191  mXCoord0=mpZScatt->GetXCoord();
192  mYCoord0=mpZScatt->GetYCoord();
193  mZCoord0=mpZScatt->GetZCoord();
194  if(mAtomWeight.numElements() != mXCoord0.numElements())
195  {
196  mAtomWeight.resize(mXCoord0.numElements());
197  mAtomWeight=1;
198  }
199 }
200 void ZMoveMinimizer::SetZAtomWeight(const CrystVector_REAL weight) {mAtomWeight=weight;}
201 void ZMoveMinimizer::MinimizeChange(long nbTrial)
202 {
203  if(mAtomWeight.max()<1e-3) return;
204  mOptimObj.Optimize(nbTrial,true);
205 }
206 
207 //######################################################################
208 //
209 // ZScatterer
210 //
211 //
212 //######################################################################
213 
214 ZScatterer::ZScatterer(const string &name,Crystal &cryst,
215  const REAL x,const REAL y,const REAL z,
216  const REAL phi,const REAL chi, const REAL psi):
217 mScattCompList(0),mNbAtom(0),mNbDummyAtom(0),
218 mPhi(0),mChi(0),mPsi(0),
219 mZAtomRegistry("List of ZAtoms"),
220 mCenterAtomIndex(0),
221 mPhiChiPsiMatrix(3,3),
222 mUseGlobalScattPow(false),mpGlobalScattPow(0),
223 mpZMoveMinimizer(0)
224 {
225  VFN_DEBUG_MESSAGE("ZScatterer::ZScatterer():("<<mName<<")",5)
226  mName=name;
227  mXYZ(0)=x;
228  mXYZ(1)=y;
229  mXYZ(2)=z;
230  mPhi=phi;
231  mChi=chi;
232  mPsi=psi;
233  this->SetCrystal(cryst);
234  this->InitRefParList();
235  VFN_DEBUG_MESSAGE("ZScatterer::ZScatterer():("<<mName<<")",5)
236 }
237 
238 ZScatterer::ZScatterer(const ZScatterer &old):
239 Scatterer(old),m3DDisplayIndex(old.m3DDisplayIndex),
240 mNbAtom(0),mNbDummyAtom(0),
241 mPhi(old.mPhi),mChi(old.mChi),mPsi(old.mPsi),
242 mCenterAtomIndex(old.mCenterAtomIndex),
243 mPhiChiPsiMatrix(old.mPhiChiPsiMatrix),
244 mUseGlobalScattPow(false),mpGlobalScattPow(0),
245 mpZMoveMinimizer(0)
246 {
247  VFN_DEBUG_ENTRY("ZScatterer::ZScatterer(&old):("<<mName<<")",10)
248 
249  mName=old.GetName();
250  mXYZ(0)=old.GetX();
251  mXYZ(1)=old.GetY();
252  mXYZ(2)=old.GetZ();
253  mPhi=old.GetPhi();
254  mChi=old.GetChi();
255  mPsi=old.GetPsi();
256  this->SetCrystal(*(old.mpCryst));
257 
258  this->InitRefParList();
259 
260  VFN_DEBUG_MESSAGE("ZScatterer::ZScatterer(&old):Copying atoms",10)
261  for(long i=0; i<old.GetZAtomRegistry().GetNb();i++)
262  this->AddAtom(old.GetZAtomRegistry().GetObj(i).GetName(),
263  old.GetZAtomRegistry().GetObj(i).GetScatteringPower(),
264  old.GetZAtomRegistry().GetObj(i).GetZBondAtom(),
265  old.GetZAtomRegistry().GetObj(i).GetZBondLength(),
266  old.GetZAtomRegistry().GetObj(i).GetZAngleAtom(),
267  old.GetZAtomRegistry().GetObj(i).GetZAngle(),
268  old.GetZAtomRegistry().GetObj(i).GetZDihedralAngleAtom(),
269  old.GetZAtomRegistry().GetObj(i).GetZDihedralAngle(),
270  old.GetZAtomRegistry().GetObj(i).GetOccupancy());
271 
272  this->SetUseGlobalScatteringPower(old.mUseGlobalScattPow);
273 
274  // Copy parameters attributes (limits, etc...)
275  VFN_DEBUG_MESSAGE("ZScatterer::ZScatterer(&old):Copying param attributes",10)
276  this->GetPar(mXYZ.data()). CopyAttributes(old.GetPar(old.mXYZ.data()));
277  this->GetPar(mXYZ.data()+1).CopyAttributes(old.GetPar(old.mXYZ.data()+1));
278  this->GetPar(mXYZ.data()+2).CopyAttributes(old.GetPar(old.mXYZ.data()+2));
279  this->GetPar(&mOccupancy). CopyAttributes(old.GetPar(&(old.mOccupancy)));
280  this->GetPar(&mPhi). CopyAttributes(old.GetPar(&(old.mPhi)));
281  this->GetPar(&mChi). CopyAttributes(old.GetPar(&(old.mChi)));
282  this->GetPar(&mPsi). CopyAttributes(old.GetPar(&(old.mPsi)));
283  VFN_DEBUG_MESSAGE("ZScatterer::ZScatterer(&old):Copying atoms param attributes",10)
284  for(long i=0; i<old.GetZAtomRegistry().GetNb();i++)
285  {
286  this->GetPar(&(this->GetZAtomRegistry().GetObj(i).mBondLength)).
287  CopyAttributes(old.GetPar(&(old.GetZAtomRegistry().GetObj(i).mBondLength)));
288  this->GetPar(&(this->GetZAtomRegistry().GetObj(i).mAngle)).
289  CopyAttributes(old.GetPar(&(old.GetZAtomRegistry().GetObj(i).mAngle)));
290  this->GetPar(&(this->GetZAtomRegistry().GetObj(i).mDihed)).
291  CopyAttributes(old.GetPar(&(old.GetZAtomRegistry().GetObj(i).mDihed)));
292  this->GetPar(&(this->GetZAtomRegistry().GetObj(i).mOccupancy)).
293  CopyAttributes(old.GetPar(&(old.GetZAtomRegistry().GetObj(i).mOccupancy)));
294  }
295 
296  VFN_DEBUG_EXIT("ZScatterer::ZScatterer(&old):("<<mName<<")",10)
297 }
298 
299 ZScatterer::~ZScatterer()
300 {
301  VFN_DEBUG_MESSAGE("ZScatterer::~ZScatterer():("<<mName<<")",5)
302  if(0 != mpGlobalScattPow) delete mpGlobalScattPow;
303  mZAtomRegistry.DeleteAll();
304 }
305 
306 ZScatterer* ZScatterer::CreateCopy() const
307 {
308  VFN_DEBUG_MESSAGE("ZScatterer::CreateCopy()"<<mName<<")",5)
309  return new ZScatterer(*this);
310 }
311 const string& ZScatterer::GetClassName() const
312 {
313  const static string className="ZScatterer";
314  return className;
315 }
316 
317 void ZScatterer::AddAtom(const string &name,const ScatteringPower *pow,
318  const long atomBond, const REAL bondLength,
319  const long atomAngle, const REAL bondAngle,
320  const long atomDihedral, const REAL dihedralAngle,
321  const REAL popu)
322 {
323  VFN_DEBUG_MESSAGE("ZScatterer::AddAtom():"<<name<<"):"<<atomBond<<" / "<<atomAngle<<" / "
324  <<atomDihedral<<" / "<<bondLength<<","<<bondAngle<<","<<dihedralAngle,10)
325  ZAtom *zatom =new ZAtom(*this,pow,
326  atomBond,bondLength,
327  atomAngle,bondAngle,
328  atomDihedral,dihedralAngle,
329  popu,name);
330 
331  if(true==mUseGlobalScattPow)
332  {
333  throw ObjCrystException("ZScatterer::AddAtom(() Cannot add an atom ! \
334 You are using the Global ScatteringPower approximation !!");
335 
336  }
337  bool usedBond=true,usedAngle=true,usedDihed=true;
338  if(mZAtomRegistry.GetNb()<1) usedBond=false;
339  if(mZAtomRegistry.GetNb()<2) usedAngle=false;
340  if(mZAtomRegistry.GetNb()<3) usedDihed=false;
341  char buf [20];
342  {
343  sprintf(buf,"%d-%d",(int)mNbAtom,(int)atomBond);
344  RefinablePar tmp("Length"+(string)buf,&(zatom->mBondLength),
345  1.,5.,
346 // bondLength*.9,bondLength*1.1,
347  gpRefParTypeScattConformBondLength,
348  REFPAR_DERIV_STEP_ABSOLUTE,true,false,usedBond,false,1.);
349  tmp.AssignClock(mClockScatterer);
350  this->AddPar(tmp);
351  }
352  {
353  sprintf(buf,"%d-%d-%d",(int)mNbAtom,(int)atomBond,(int)atomAngle);
354  RefinablePar tmp("Angle"+(string)buf,&(zatom->mAngle),
355  0,2*M_PI,
356 // zatom->mAngle-.2,zatom->mAngle+.2,
357  gpRefParTypeScattConformBondAngle,
358  REFPAR_DERIV_STEP_ABSOLUTE,false,false,usedAngle,true,RAD2DEG,2*M_PI);
359  tmp.AssignClock(mClockScatterer);
360  this->AddPar(tmp);
361  }
362  {
363  sprintf(buf,"%d-%d-%d-%d",(int)mNbAtom,(int)atomBond,(int)atomAngle,(int)atomDihedral);
364  RefinablePar tmp("Dihed"+(string)buf,&(zatom->mDihed),
365  0,2*M_PI,
366 // zatom->mDihed-.2,zatom->mDihed+.2,
367  gpRefParTypeScattConformDihedAngle,
368  REFPAR_DERIV_STEP_ABSOLUTE,false,false,usedDihed,true,RAD2DEG,2*M_PI);
369  tmp.AssignClock(mClockScatterer);
370  this->AddPar(tmp);
371  }
372  {//fixed by default
373  sprintf(buf,"%d",(int)mNbAtom);
374  RefinablePar tmp("Occupancy"+(string)buf,
375  &(zatom->mOccupancy),0,1,
376  gpRefParTypeScattOccup,
377  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,1.,1.);
378  tmp.AssignClock(mClockScatterer);
379  this->AddPar(tmp);
380  }
381 
382  //:NOTE: we *must* add it in the registry after declaring the parameters,
383  // since it triggers the creation of the WXZAtom, which requires the parameters...
384  VFN_DEBUG_MESSAGE("ZScatterer::AddAtom():Registering...",2)
385  mZAtomRegistry.Register(*zatom);
386 
387  // Add an entry for this atom in the list of all components. The actual values are
388  // written in Update().No entry for a dummy atom
389  VFN_DEBUG_MESSAGE("ZScatterer::AddAtom():Adding to the ScattCompList...",2)
390  if(0!=pow)
391  {
392  ++mScattCompList;
393  mScattCompList(mNbAtom-mNbDummyAtom).mpScattPow=pow;
394  mComponentIndex.resizeAndPreserve(mNbAtom-mNbDummyAtom+1);
395  mComponentIndex(mNbAtom-mNbDummyAtom)=mNbAtom;
396  }
397  else mNbDummyAtom++;
398 
399  //Finish
400  mNbAtom++;
401  mClockScatterer.Click();
402  VFN_DEBUG_MESSAGE("ZScatterer::AddAtom():End",3)
403 }
404 
405 int ZScatterer::GetNbComponent() const
406 {
407  if(true==mUseGlobalScattPow) return 1;
408  return mNbAtom-mNbDummyAtom;
409 }
410 const ScatteringComponentList& ZScatterer::GetScatteringComponentList() const
411 {
412  this->UpdateScattCompList();
413 
414  return mScattCompList;
415 }
416 
417 string ZScatterer::GetComponentName(const int i) const
418 {
419  return mZAtomRegistry.GetObj(mComponentIndex(i)).GetName();
420 }
421 
422 void ZScatterer::Print() const
423 {
424  VFN_DEBUG_MESSAGE("ZScatterer::Print()",5)
425  //:TODO:
426  //this->UpdateScattCompList();
427  //for(int i=0;i<this->mNbAtom;i++) ??;
428 }
429 
430 REAL ZScatterer::GetPhi() const {return mPhi;}
431 REAL ZScatterer::GetChi() const {return mChi;}
432 REAL ZScatterer::GetPsi() const {return mPsi;}
433 
434 void ZScatterer::SetPhi(const REAL x) { mClockScatterer.Click();mPhi=x;}
435 void ZScatterer::SetChi(const REAL y) { mClockScatterer.Click();mChi=y;}
436 void ZScatterer::SetPsi(const REAL z) { mClockScatterer.Click();mPsi=z;}
437 
438 REAL ZScatterer::GetZAtomX(const int i)const{this->UpdateCoordinates(); return mXCoord(i);}
439 REAL ZScatterer::GetZAtomY(const int i)const{this->UpdateCoordinates(); return mYCoord(i);}
440 REAL ZScatterer::GetZAtomZ(const int i)const{this->UpdateCoordinates(); return mZCoord(i);}
441 long ZScatterer::GetZBondAtom(const int i)const
442 {return mZAtomRegistry.GetObj(i).GetZBondAtom();}
443 
444 long ZScatterer::GetZAngleAtom(const int i)const
445 {return mZAtomRegistry.GetObj(i).GetZAngleAtom();}
446 
447 long ZScatterer::GetZDihedralAngleAtom(const int i)const
448 {return mZAtomRegistry.GetObj(i).GetZDihedralAngleAtom();}
449 
450 REAL ZScatterer::GetZBondLength(const int i)const
451 {return mZAtomRegistry.GetObj(i).GetZBondLength();}
452 REAL ZScatterer::GetZAngle(const int i)const
453 {return mZAtomRegistry.GetObj(i).GetZAngle();}
454 REAL ZScatterer::GetZDihedralAngle(const int i)const
455 {return mZAtomRegistry.GetObj(i).GetZDihedralAngle();}
456 
457 void ZScatterer::SetZBondLength(const int i,const REAL a)
458 {mClockScatterer.Click();mZAtomRegistry.GetObj(i).SetZBondLength(a);}
459 
460 void ZScatterer::SetZAngle(const int i,const REAL a)
461 {mClockScatterer.Click();mZAtomRegistry.GetObj(i).SetZAngle(a);}
462 
463 void ZScatterer::SetZDihedralAngle(const int i,const REAL a)
464  {mClockScatterer.Click();mZAtomRegistry.GetObj(i).SetZDihedralAngle(a);}
465 
466 const ObjRegistry<ZAtom>& ZScatterer::GetZAtomRegistry()const
467 {return mZAtomRegistry;}
468 
469 ostream& ZScatterer::POVRayDescription(ostream &os,
470  const CrystalPOVRayOptions &options)const
471 {
472  VFN_DEBUG_MESSAGE("ZScatterer::POVRayDescription()",5)
473  //throw ObjCrystException("ZScatterer::POVRayDescription() not implemented! "+this->GetName());
474  //:TODO:
475  this->UpdateScattCompList();
476 
477  //for(long i=0;i<mNbAtom;i++) mpAtom[i]->POVRayDescription(os,onlyIndependentAtoms);
478  os << "// Description of ZScatterer :" << this->GetName() <<endl;
479 
480 #if 0
481  if(true==mUseGlobalScattPow)
482  {
483  mpAtom[mCenterAtomIndex]->POVRayDescription(os,onlyIndependentAtoms);
484  return os;
485  }
486  //Declare colours
487  for(int i=0;i<mNbAtom;i++)
488  {
489  if(mpAtom[i]->IsDummy()) continue;
490  os <<" #declare colour_"<<mpAtom[i]->GetName()<<"="<<mpAtom[i]->Colour()<<";"<<endl;
491  }
492 
493  if(true==onlyIndependentAtoms)
494  {
495  CrystVector_REAL x(mNbAtom),y(mNbAtom),z(mNbAtom);
496  for(int i=0;i<mNbAtom;i++)
497  {
498  x(i)=mpAtom[i]->GetX();
499  y(i)=mpAtom[i]->GetY();
500  z(i)=mpAtom[i]->GetZ();
501  this->GetCrystal().FractionalToOrthonormalCoords(x(i),y(i),z(i));
502  }
503  for(int i=0;i<mNbAtom;i++)
504  {
505  if(mpAtom[i]->IsDummy())
506  {
507  os << " // Skipping Dummy Atom :" << mpAtom[i]->GetName() <<endl<<endl;
508  continue;
509  }
510  os << " // Atom :" << mpAtom[i]->GetName() <<endl;
511  os << " sphere " << endl;
512  os << " { <"<<x(i)<<","<<y(i)<<","<<z(i)<< ">,"
513  << mpAtom[i]->GetRadius()/3<<endl;
514  os << " finish { ambient 0.2 diffuse 0.8 phong 1}" <<endl;
515  os << " pigment { colour colour_"<< mpAtom[i]->GetName() <<" }" << endl;
516  os << " }" <<endl;
517  //Draw the bond for this Atom,if it's not linked to a dummy
518  int bond=ZBondAtom(i);
519  if((mpAtom[bond]->IsDummy()) || (i==0)) continue;
520  os << " cylinder"<<endl;
521  os << " { <"<<x(i)<<","<<y(i)<<","<<z(i)<< ">,"<<endl;
522  os << " <"<<x(bond)<<","<<y(bond)<<","<<z(bond)<< ">,"<<endl;
523  os << " 0.1"<<endl;
524  os << " pigment { colour Gray }"<<endl;
525  os << " }"<<endl;
526  }
527  }
528  else
529  {
530  vector<CrystMatrix_REAL> xyzCoords;
531  for(int i=0;i<mNbAtom;i++)
532  vXYZCoords.push_back(this->GetCrystal().GetSpaceGroup().GetAllSymmetrics(mpAtom[i]->GetX(),
533  mpAtom[i]->GetY(),
534  mpAtom[i]->GetZ(),
535  false,false,false));
536  const int nbSymmetrics=(vXYZCoords[0])->rows();
537  int symNum=0;
538  CrystMatrix_int translate(27,3);
539  translate= -1,-1,-1,
540  -1,-1, 0,
541  -1,-1, 1,
542  -1, 0,-1,
543  -1, 0, 0,
544  -1, 0, 1,
545  -1, 1,-1,
546  -1, 1, 0,
547  -1, 1, 1,
548  0,-1,-1,
549  0,-1, 0,
550  0,-1, 1,
551  0, 0,-1,
552  0, 0, 0,
553  0, 0, 1,
554  0, 1,-1,
555  0, 1, 0,
556  0, 1, 1,
557  1,-1,-1,
558  1,-1, 0,
559  1,-1, 1,
560  1, 0,-1,
561  1, 0, 0,
562  1, 0, 1,
563  1, 1,-1,
564  1, 1, 0,
565  1, 1, 1;
566  REAL dx,dy,dz;
567  CrystVector_REAL x(mNbAtom),y(mNbAtom),z(mNbAtom);
568  CrystVector_REAL xSave,ySave,zSave;
569  for(int i=0;i<nbSymmetrics;i++)
570  {
571  for(int j=0;j<mNbAtom;j++)
572  {
573  x(j)=vXYZCoords[j](i,0);
574  y(j)=vXYZCoords[j](i,1);
575  z(j)=vXYZCoords[j](i,2);
576  }
577  //Bring back central atom in unit cell; move peripheral atoms with the same amount
578  dx=x(0);
579  dy=y(0);
580  dz=z(0);
581  x(0) = fmod((float) x(0),(float)1); if(x(0)<0) x(0)+=1.;
582  y(0) = fmod((float) y(0),(float)1); if(y(0)<0) y(0)+=1.;
583  z(0) = fmod((float) z(0),(float)1); if(z(0)<0) z(0)+=1.;
584  dx = x(0)-dx;
585  dy = y(0)-dy;
586  dz = z(0)-dz;
587  for(int j=1;j<mNbAtom;j++)
588  {
589  x(j) += dx;
590  y(j) += dy;
591  z(j) += dz;
592  }
593  //Generate also translated atoms near the unit cell
594  const REAL limit =0.1;
595  for(int j=0;j<translate.rows();j++)
596  {
597  xSave=x;
598  ySave=y;
599  zSave=z;
600  x += translate(j,0);
601  y += translate(j,1);
602  z += translate(j,2);
603  if( (x(0)>(-limit)) && (x(0)<(1+limit))
604  &&(y(0)>(-limit)) && (y(0)<(1+limit))
605  &&(z(0)>(-limit)) && (z(0)<(1+limit)))
606  {
607  for(int k=0;k<mNbAtom;k++)
608  this->GetCrystal().FractionalToOrthonormalCoords(x(k),y(k),z(k));
609  os << " // Symmetric&Translated #" << symNum++ <<endl;
610  //:NOTE: The code below is the same as without symmetrics
611  for(int i=0;i<mNbAtom;i++)
612  {
613  if(mpAtom[i]->IsDummy())
614  {
615  os << " // Skipping Dummy Atom :" << mpAtom[i]->GetName() <<endl<<endl;
616  continue;
617  }
618  os << " // Atom :" << mpAtom[i]->GetName() <<endl;
619  os << " sphere " << endl;
620  os << " { <"<<x(i)<<","<<y(i)<<","<<z(i)<< ">,"
621  << mpAtom[i]->GetRadius()/3<<endl;
622  os << " finish { ambient 0.2 diffuse 0.8 phong 1}" <<endl;
623  os << " pigment { colour colour_"<< mpAtom[i]->GetName() <<" }" << endl;
624  os << " }" <<endl;
625  //Draw the bond for this Atom,if it's not linked to a dummy
626  int bond=ZBondAtom(i);
627  if((mpAtom[bond]->IsDummy()) || (i==0)) continue;
628  os << " cylinder"<<endl;
629  os << " { <"<<x(i)<<","<<y(i)<<","<<z(i)<< ">,"<<endl;
630  os << " <"<<x(bond)<<","<<y(bond)<<","<<z(bond)<< ">,"<<endl;
631  os << " 0.1"<<endl;
632  os << " pigment { colour Gray }"<<endl;
633  os << " }"<<endl;
634  }
635  }//if in limits
636  x=xSave;
637  y=ySave;
638  z=zSave;
639  }//for translation
640  }//for symmetrics
641  }//else
642 #endif
643  return os;
644 }
645 
646 #ifdef OBJCRYST_GL
647 void ZScatterer::GLInitDisplayList(const bool onlyIndependentAtoms,
648  const REAL xMin,const REAL xMax,
649  const REAL yMin,const REAL yMax,
650  const REAL zMin,const REAL zMax,
651  const bool displayEnantiomer,
652  const bool displayNames,
653  const bool hideHydrogens,
654  const REAL fadeDistance,
655  const bool fullMoleculeInLimits)const
656 {
657  VFN_DEBUG_ENTRY("ZScatterer::GLInitDisplayList()",4)
658  if(mZAtomRegistry.GetNb()==0)
659  {
660  VFN_DEBUG_EXIT("ZScatterer::GLInitDisplayList():No ZAtom to display !",4)
661  return;
662  }
663  REAL en=1;
664  if(displayEnantiomer==true) en=-1;
665  this->UpdateScattCompList();
666  if(true==mUseGlobalScattPow)
667  {
668  //mpAtom[mCenterAtomIndex]->GLInitDisplayList(onlyIndependentAtoms);
669  return;
670  }
671 
672  GLfloat colour_bond[]= { 0.5, .5, .5, 1.0 };
673  GLfloat colour_side[]= { 0.0, .0, .0, 1.0 };
674  const GLfloat colour0[] = {0.0f, 0.0f, 0.0f, 0.0f};
675 
676  GLUquadricObj* pQuadric = gluNewQuadric();
677 
678  if(true==onlyIndependentAtoms)
679  {
680  //cout << m3DDisplayIndex <<endl;
681  CrystVector_REAL x,y,z;
682  x=mXCoord;
683  y=mYCoord;
684  z=mZCoord;
685  if(m3DDisplayIndex.numElements()>0)
686  {
687  for(long k=0;k<m3DDisplayIndex.rows();k++)
688  {
689  switch(m3DDisplayIndex(k,0))
690  {
691  case 0:break;
692  case 1: //Draw a sphere
693  {
694  const int n1=m3DDisplayIndex(k,1);
695  if(0==mZAtomRegistry.GetObj(n1).GetScatteringPower())break;
696  if(hideHydrogens && (mZAtomRegistry.GetObj(n1).GetScatteringPower()->GetForwardScatteringFactor(RAD_XRAY)<1.5)) continue;
697  glMaterialfv (GL_FRONT, GL_AMBIENT_AND_DIFFUSE,
698  mZAtomRegistry.GetObj(n1).GetScatteringPower()->GetColourRGB());
699  glPushMatrix();
700  glTranslatef(x(n1)*en, y(n1), z(n1));
701  gluSphere(pQuadric,mZAtomRegistry.GetObj(n1).GetScatteringPower()
702  ->GetRadius()/3.,10,10);
703  glPopMatrix();
704  }
705  case 2: // Draw a bond
706  {
707  long n1,n2;
708  n1=m3DDisplayIndex(k,1);
709  n2=m3DDisplayIndex(k,2);
710  if(hideHydrogens && (mZAtomRegistry.GetObj(n1).GetScatteringPower()->GetForwardScatteringFactor(RAD_XRAY)<1.5)) continue;
711  if(hideHydrogens && (mZAtomRegistry.GetObj(n2).GetScatteringPower()->GetForwardScatteringFactor(RAD_XRAY)<1.5)) continue;
712  glPushMatrix();
713  glTranslatef(x(n1)*en, y(n1), z(n1));
714  glMaterialfv (GL_FRONT, GL_AMBIENT_AND_DIFFUSE,colour_bond);
715  GLUquadricObj *quadobj = gluNewQuadric();
716  glColor3f(1.0f,1.0f,1.0f);
717  const REAL height= sqrt( (x(n2)-x(n1))*(x(n2)-x(n1))
718  +(y(n2)-y(n1))*(y(n2)-y(n1))
719  +(z(n2)-z(n1))*(z(n2)-z(n1)));
720  glRotatef(180,(x(n2)-x(n1))*en,y(n2)-y(n1),z(n2)-z(n1)+height);// ?!?!?!
721  gluCylinder(quadobj,.1,.1,height,10,1 );
722  gluDeleteQuadric(quadobj);
723  glPopMatrix();
724 
725  }
726  case 3: // Draw a triangular face
727  {
728  long n1,n2,n3;
729  REAL x1,y1,z1,x2,y2,z2,xn,yn,zn,xc,yc,zc;
730  n1=m3DDisplayIndex(k,1);
731  n2=m3DDisplayIndex(k,2);
732  n3=m3DDisplayIndex(k,3);
733  x1=x(n1)-x(n2);
734  y1=y(n1)-y(n2);
735  z1=z(n1)-z(n2);
736 
737  x2=x(n1)-x(n3);
738  y2=y(n1)-y(n3);
739  z2=z(n1)-z(n3);
740 
741  xn=y1*z2-z1*y2;
742  yn=z1*x2-x1*z2;
743  zn=x1*y2-y1*x2;
744 
745  xc=(x(n1)+x(n2)+x(n3))/3.-x(mCenterAtomIndex);
746  yc=(y(n1)+y(n2)+y(n3))/3.-y(mCenterAtomIndex);
747  zc=(z(n1)+z(n2)+z(n3))/3.-z(mCenterAtomIndex);
748 
749  glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,
750  mZAtomRegistry.GetObj(0).GetScatteringPower()->GetColourRGB());
751  //glColor3f(1.0f,1.0f,1.0f); // White
752  glBegin(GL_TRIANGLES); // Bottom
753  if((xn*xc+yn*yc+zn*zc)>0) glNormal3f(xn*en, yn, zn);
754  else glNormal3f(-xn*en, -yn, -zn);
755  glVertex3f(x(n1)*en,y(n1),z(n1));
756  glVertex3f(x(n2)*en,y(n2),z(n2));
757  glVertex3f(x(n3)*en,y(n3),z(n3));
758  glEnd();
759  glMaterialfv (GL_FRONT, GL_AMBIENT_AND_DIFFUSE,colour_side);
760  glBegin(GL_LINE_LOOP);
761  glVertex3f(x(n1)*en,y(n1),z(n1));
762  glVertex3f(x(n2)*en,y(n2),z(n2));
763  glVertex3f(x(n3)*en,y(n3),z(n3));
764  glEnd();
765  }
766  case 4: // Draw a quadric face
767  {
768  long n1,n2,n3,n4;
769  REAL x1,y1,z1,x2,y2,z2,xn,yn,zn,xc,yc,zc;
770  n1=m3DDisplayIndex(k,1);
771  n2=m3DDisplayIndex(k,2);
772  n3=m3DDisplayIndex(k,3);
773  n4=m3DDisplayIndex(k,3);
774  x1=x(n1)-x(n2);
775  y1=y(n1)-y(n2);
776  z1=z(n1)-z(n2);
777 
778  x2=x(n1)-x(n3);
779  y2=y(n1)-y(n3);
780  z2=z(n1)-z(n3);
781 
782  xn=y1*z2-z1*y2;
783  yn=z1*x2-x1*z2;
784  zn=x1*y2-y1*x2;
785 
786  xc=(x(n1)+x(n2)+x(n3)+x(n4))/4.-x(mCenterAtomIndex);
787  yc=(y(n1)+y(n2)+y(n3)+y(n4))/4.-y(mCenterAtomIndex);
788  zc=(z(n1)+z(n2)+z(n3)+z(n4))/4.-z(mCenterAtomIndex);
789 
790  glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,
791  mZAtomRegistry.GetObj(0).GetScatteringPower()->GetColourRGB());
792  //glColor3f(1.0f,1.0f,1.0f); // White
793  glBegin(GL_TRIANGLES); // Bottom
794  if((xn*xc+yn*yc+zn*zc)>0) glNormal3f(xn*en, yn, zn);
795  else glNormal3f(-xn*en, -yn, -zn);
796  glVertex3f(x(n1)*en,y(n1),z(n1));
797  glVertex3f(x(n2)*en,y(n2),z(n2));
798  glVertex3f(x(n3)*en,y(n3),z(n3));
799  glVertex3f(x(n4)*en,y(n4),z(n4));
800  glEnd();
801  glMaterialfv (GL_FRONT, GL_AMBIENT_AND_DIFFUSE,colour_side);
802  glBegin(GL_LINE_LOOP);
803  glVertex3f(x(n1)*en,y(n1),z(n1));
804  glVertex3f(x(n2)*en,y(n2),z(n2));
805  glVertex3f(x(n3)*en,y(n3),z(n3));
806  glVertex3f(x(n4)*en,y(n4),z(n4));
807  glEnd();
808  }
809  }
810  }
811  }
812  else
813  {
814  for(int k=0;k<mNbAtom;k++)
815  {
816  if(0==mZAtomRegistry.GetObj(k).GetScatteringPower())continue;
817  if(hideHydrogens && (mZAtomRegistry.GetObj(k).GetScatteringPower()->GetForwardScatteringFactor(RAD_XRAY)<1.5)) continue;
818  const float r=mZAtomRegistry.GetObj(k).GetScatteringPower()->GetColourRGB()[0];
819  const float g=mZAtomRegistry.GetObj(k).GetScatteringPower()->GetColourRGB()[1];
820  const float b=mZAtomRegistry.GetObj(k).GetScatteringPower()->GetColourRGB()[2];
821  glPushMatrix();
822  glTranslatef(x(k)*en, y(k), z(k));
823  if(displayNames)
824  {
825  GLfloat colourChar [] = {1.0, 1.0, 1.0, 1.0};
826  if((r>0.8)&&(g>0.8)&&(b>0.8))
827  {
828  colourChar[0] = 0.5;
829  colourChar[1] = 0.5;
830  colourChar[2] = 0.5;
831  }
832  glMaterialfv(GL_FRONT, GL_AMBIENT, colour0);
833  glMaterialfv(GL_FRONT, GL_DIFFUSE, colour0);
834  glMaterialfv(GL_FRONT, GL_SPECULAR, colour0);
835  glMaterialfv(GL_FRONT, GL_EMISSION, colourChar);
836  glMaterialfv(GL_FRONT, GL_SHININESS,colour0);
837  glRasterPos3f(0.0f, 0.0f, 0.0f);
838  crystGLPrint(mZAtomRegistry.GetObj(k).GetName());
839  }
840  else
841  {
842  const GLfloat colourAtom [] = {r, g, b, 1.0};
843  glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE,colourAtom);
844  glMaterialfv(GL_FRONT, GL_SPECULAR, colour0);
845  glMaterialfv(GL_FRONT, GL_EMISSION, colour0);
846  glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
847  glPolygonMode(GL_FRONT, GL_FILL);
848  gluSphere(pQuadric,
849  mZAtomRegistry.GetObj(k).GetScatteringPower()->GetRadius()/3.,10,10);
850  //Draw the bond for this Atom,if it's not linked to a dummy
851  int bond=this->GetZBondAtom(k);
852  if((0!=mZAtomRegistry.GetObj(bond).GetScatteringPower()) && (k>0))
853  {
854  glMaterialfv(GL_FRONT,GL_AMBIENT_AND_DIFFUSE,colour_bond);
855  GLUquadricObj *quadobj = gluNewQuadric();
856  glColor3f(1.0f,1.0f,1.0f);
857  const REAL height= sqrt( (x(bond)-x(k))*(x(bond)-x(k))
858  +(y(bond)-y(k))*(y(bond)-y(k))
859  +(z(bond)-z(k))*(z(bond)-z(k)));
860  glRotatef(180,(x(bond)-x(k))*en,y(bond)-y(k),z(bond)-z(k)+height);// !!!
861  gluCylinder(quadobj,.1,.1,height,10,1 );
862  gluDeleteQuadric(quadobj);
863  }
864  }
865  glPopMatrix();
866  }
867  }//Use triangle faces ?
868  }//Only independent atoms ?
869  else
870  {
871  VFN_DEBUG_ENTRY("ZScatterer::GLInitDisplayList():Show all symmetrics",3)
872  vector<CrystMatrix_REAL> vXYZCoords;
873  {
874  REAL x0,y0,z0;
875  for(int i=0;i<mNbAtom;i++)
876  {//We also generate the positions of dummy atoms.. This may be needed...
877  x0=mXCoord(i);
878  y0=mYCoord(i);
879  z0=mZCoord(i);
880  this->GetCrystal().OrthonormalToFractionalCoords(x0,y0,z0);
881  vXYZCoords.push_back(this->GetCrystal().GetSpaceGroup().
882  GetAllSymmetrics(x0,y0,z0,false,false,false));
883  }
884  }
885  CrystMatrix_int translate(27,3);
886  translate= -1,-1,-1,
887  -1,-1, 0,
888  -1,-1, 1,
889  -1, 0,-1,
890  -1, 0, 0,
891  -1, 0, 1,
892  -1, 1,-1,
893  -1, 1, 0,
894  -1, 1, 1,
895  0,-1,-1,
896  0,-1, 0,
897  0,-1, 1,
898  0, 0,-1,
899  0, 0, 0,
900  0, 0, 1,
901  0, 1,-1,
902  0, 1, 0,
903  0, 1, 1,
904  1,-1,-1,
905  1,-1, 0,
906  1,-1, 1,
907  1, 0,-1,
908  1, 0, 0,
909  1, 0, 1,
910  1, 1,-1,
911  1, 1, 0,
912  1, 1, 1;
913  REAL dx,dy,dz;
914  CrystVector_REAL x(mNbAtom),y(mNbAtom),z(mNbAtom);
915  CrystVector_REAL xSave,ySave,zSave;
916  const int nbSymmetrics=vXYZCoords[0].rows();
917  for(int i=0;i<nbSymmetrics;i++)
918  {
919  VFN_DEBUG_ENTRY("ZScatterer::GLInitDisplayList():Symmetric#"<<i,3)
920  for(int j=0;j<mNbAtom;j++)
921  {
922  x(j)=vXYZCoords[j](i,0);
923  y(j)=vXYZCoords[j](i,1);
924  z(j)=vXYZCoords[j](i,2);
925  }
926  //Bring back central atom in unit cell; move peripheral atoms with the same amount
927  dx=x(0);
928  dy=y(0);
929  dz=z(0);
930  x(0) = fmod((float) x(0),(float)1); if(x(0)<0) x(0)+=1.;
931  y(0) = fmod((float) y(0),(float)1); if(y(0)<0) y(0)+=1.;
932  z(0) = fmod((float) z(0),(float)1); if(z(0)<0) z(0)+=1.;
933  dx = x(0)-dx;
934  dy = y(0)-dy;
935  dz = z(0)-dz;
936  for(int j=1;j<mNbAtom;j++)
937  {
938  x(j) += dx;
939  y(j) += dy;
940  z(j) += dz;
941  }
942  //Generate also translated atoms near the unit cell
943  xSave=x;
944  ySave=y;
945  zSave=z;
946  for(int j=0;j<translate.rows();j++)
947  {
948  x += translate(j,0);
949  y += translate(j,1);
950  z += translate(j,2);
951  if( (x(0)>xMin) && (x(0)<xMax)
952  &&(y(0)>yMin) && (y(0)<yMax)
953  &&(z(0)>zMin) && (z(0)<zMax))
954  {
955  for(int k=0;k<mNbAtom;k++)
956  this->GetCrystal().FractionalToOrthonormalCoords(x(k),y(k),z(k));
957  //:NOTE: The code below is the same as without symmetrics
958  if(m3DDisplayIndex.numElements()>0)
959  {
960  long n1,n2,n3;
961  REAL x1,y1,z1,x2,y2,z2,xn,yn,zn,xc,yc,zc;
962  for(long k=0;k<m3DDisplayIndex.rows();k++)
963  {
964  n1=m3DDisplayIndex(k,0);
965  n2=m3DDisplayIndex(k,1);
966  n3=m3DDisplayIndex(k,2);
967  x1=x(n1)-x(n2);
968  y1=y(n1)-y(n2);
969  z1=z(n1)-z(n2);
970 
971  x2=x(n1)-x(n3);
972  y2=y(n1)-y(n3);
973  z2=z(n1)-z(n3);
974 
975  xn=y1*z2-z1*y2;
976  yn=z1*x2-x1*z2;
977  zn=x1*y2-y1*x2;
978 
979  xc=(x(n1)+x(n2)+x(n3))/3.-x(mCenterAtomIndex);
980  yc=(y(n1)+y(n2)+y(n3))/3.-y(mCenterAtomIndex);
981  zc=(z(n1)+z(n2)+z(n3))/3.-z(mCenterAtomIndex);
982 
983  glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,
984  mZAtomRegistry.GetObj(0).GetScatteringPower()->GetColourRGB());
985  glBegin(GL_TRIANGLES);
986  if((xn*xc+yn*yc+zn*zc)>0) glNormal3f( xn*en, yn, zn);
987  else glNormal3f(-xn*en, -yn, -zn);
988 
989  glVertex3f(x(n1)*en,y(n1),z(n1));
990  glVertex3f(x(n2)*en,y(n2),z(n2));
991  glVertex3f(x(n3)*en,y(n3),z(n3));
992  glEnd();
993  glMaterialfv (GL_FRONT, GL_AMBIENT_AND_DIFFUSE,colour_side);
994  glBegin(GL_LINE_LOOP);
995  glVertex3f(x(n1)*en,y(n1),z(n1));
996  glVertex3f(x(n2)*en,y(n2),z(n2));
997  glVertex3f(x(n3)*en,y(n3),z(n3));
998  glEnd();
999  }
1000  }
1001  else
1002  {
1003  for(int k=0;k<mNbAtom;k++)
1004  {
1005  if(0==mZAtomRegistry.GetObj(k).GetScatteringPower())continue;
1006  if(hideHydrogens && (mZAtomRegistry.GetObj(k).GetScatteringPower()->GetForwardScatteringFactor(RAD_XRAY)<1.5)) continue;
1007  const float r=mZAtomRegistry.GetObj(k).GetScatteringPower()->GetColourRGB()[0];
1008  const float g=mZAtomRegistry.GetObj(k).GetScatteringPower()->GetColourRGB()[1];
1009  const float b=mZAtomRegistry.GetObj(k).GetScatteringPower()->GetColourRGB()[2];
1010  glPushMatrix();
1011  glTranslatef(x(k)*en, y(k), z(k));
1012  if(displayNames)
1013  {
1014  GLfloat colourChar [] = {1.0, 1.0, 1.0, 1.0};
1015  if((r>0.8)&&(g>0.8)&&(b>0.8))
1016  {
1017  colourChar[0] = 0.5;
1018  colourChar[1] = 0.5;
1019  colourChar[2] = 0.5;
1020  }
1021  glMaterialfv(GL_FRONT, GL_AMBIENT, colour0);
1022  glMaterialfv(GL_FRONT, GL_DIFFUSE, colour0);
1023  glMaterialfv(GL_FRONT, GL_SPECULAR, colour0);
1024  glMaterialfv(GL_FRONT, GL_EMISSION, colourChar);
1025  glMaterialfv(GL_FRONT, GL_SHININESS,colour0);
1026  glRasterPos3f(0.0f, 0.0f, 0.0f);
1027  crystGLPrint(mZAtomRegistry.GetObj(k).GetName());
1028  }
1029  else
1030  {
1031  const GLfloat colourAtom [] = {r, g, b, 1.0};
1032  glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE,colourAtom);
1033  glMaterialfv(GL_FRONT, GL_SPECULAR, colour0);
1034  glMaterialfv(GL_FRONT, GL_EMISSION, colour0);
1035  glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
1036  glPolygonMode(GL_FRONT, GL_FILL);
1037  gluSphere(pQuadric,
1038  mZAtomRegistry.GetObj(k).GetScatteringPower()->GetRadius()/3.,10,10);
1039  glRasterPos3f(0,0,0);
1040  //Draw the bond for this Atom,if it's not linked to a dummy
1041  int bond=this->GetZBondAtom(k);
1042  if((0!=mZAtomRegistry.GetObj(bond).GetScatteringPower()) && (k>0))
1043  {
1044  glMaterialfv(GL_FRONT,GL_AMBIENT_AND_DIFFUSE,colour_bond);
1045  GLUquadricObj *quadobj = gluNewQuadric();
1046  glColor3f(1.0f,1.0f,1.0f);
1047  const REAL height= sqrt( (x(bond)-x(k))*(x(bond)-x(k))
1048  +(y(bond)-y(k))*(y(bond)-y(k))
1049  +(z(bond)-z(k))*(z(bond)-z(k)));
1050  glRotatef(180,(x(bond)-x(k))*en,y(bond)-y(k),z(bond)-z(k)+height);// !!!
1051  gluCylinder(quadobj,.1,.1,height,10,1 );
1052  gluDeleteQuadric(quadobj);
1053  }
1054  }
1055  glPopMatrix();
1056  }
1057  }//Use triangle faces ?
1058  }//if in limits
1059  x=xSave;
1060  y=ySave;
1061  z=zSave;
1062  }//for translation
1063  VFN_DEBUG_EXIT("ZScatterer::GLInitDisplayList():Symmetric#"<<i,3)
1064  }//for symmetrics
1065  VFN_DEBUG_EXIT("ZScatterer::GLInitDisplayList():Show all symmetrics",3)
1066  }//else
1067  gluDeleteQuadric(pQuadric);
1068  VFN_DEBUG_EXIT("ZScatterer::GLInitDisplayList()",4)
1069 }
1070 #endif // OBJCRYST_GL
1071 
1072 void ZScatterer::SetUseGlobalScatteringPower(const bool useIt)
1073 {
1074  VFN_DEBUG_MESSAGE("ZScatterer::SetUseGlobalScatteringPower(bool):"<<this->GetName()<<":"<<useIt,5)
1075  if(true==useIt)
1076  {
1077  if(false==mUseGlobalScattPow)
1078  {
1079  mpGlobalScattPow=new GlobalScatteringPower(*this);
1080  mUseGlobalScattPow=true;
1081  mScattCompList.Reset();
1082  ++mScattCompList;
1083  this->InitRefParList();
1084  mClockScatterer.Click();
1085  }
1086  //Just change the functions which return the component list
1087  }
1088  else
1089  {
1090  if(true==mUseGlobalScattPow)
1091  {
1092  delete mpGlobalScattPow;
1093  mpGlobalScattPow=0;
1094  mUseGlobalScattPow=false;
1095  mScattCompList.Reset();
1096  for(long i=0;i<(mNbAtom-mNbDummyAtom);i++) ++mScattCompList;
1097  this->InitRefParList();
1098  mClockScatterer.Click();
1099  }
1100  }
1101 }
1102 
1103 void ZScatterer::GetGeneGroup(const RefinableObj &obj,
1104  CrystVector_uint & groupIndex,
1105  unsigned int &first) const
1106 {
1107  // One group for all translation parameters, another for orientation.
1108  // All conformation parameters (bond, angle,torsion) are in individual groups.
1109  unsigned int posIndex=0;
1110  unsigned int orientIndex=0;
1111  VFN_DEBUG_MESSAGE("ZScatterer::GetGeneGroup()",4)
1112  for(long i=0;i<obj.GetNbPar();i++)
1113  for(long j=0;j<this->GetNbPar();j++)
1114  if(&(obj.GetPar(i)) == &(this->GetPar(j)))
1115  {
1116  if(this->GetPar(j).GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattTransl))
1117  {
1118  if(posIndex==0) posIndex=first++;
1119  groupIndex(i)=posIndex;
1120  }
1121  else
1122  if(this->GetPar(j).GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattOrient))
1123  {
1124  if(orientIndex==0) orientIndex=first++;
1125  groupIndex(i)=orientIndex;
1126  }
1127  else groupIndex(i)= first++;
1128  }
1129 }
1130 const CrystVector_REAL& ZScatterer::GetXCoord() const
1131 {
1132  this->UpdateCoordinates();
1133  return mXCoord;
1134 }
1135 const CrystVector_REAL& ZScatterer::GetYCoord() const
1136 {
1137  this->UpdateCoordinates();
1138  return mYCoord;
1139 }
1140 const CrystVector_REAL& ZScatterer::GetZCoord() const
1141 {
1142  this->UpdateCoordinates();
1143  return mZCoord;
1144 }
1145 
1146 void ZScatterer::EndOptimization()
1147 {
1148  if(mOptimizationDepth==1)
1149  {
1150  if(0!=mpZMoveMinimizer) delete mpZMoveMinimizer;
1151  mpZMoveMinimizer=0;
1152  }
1153  this->RefinableObj::EndOptimization();
1154 }
1155 
1156 std::vector<std::string> SplitString(const std::string &s)
1157 {
1158  const string delim(" ");
1159  std::vector<std::string> v;
1160  string str=s;
1161  unsigned int ct=0;
1162 
1163  int n=str.find_first_of(delim);
1164  while( n != (int) str.npos)
1165  {
1166  ct++;
1167  if(n>0)
1168  {
1169  v.push_back(str.substr(0,n));
1170  }
1171  str= str.substr(n+1);
1172  n=str.find_first_of(delim);
1173  }
1174  if(str.length() > 0) v.push_back(str);
1175 
1176  //cout<<"SplitString: "<<s<<" -> ";
1177  //for(std::vector<std::string>::const_iterator pos=v.begin();pos!=v.end();++pos) cout << *pos <<" / ";
1178  //cout<<endl;
1179 
1180  return v;
1181 }
1182 
1192 void ReadFHLine(const char*buf, const unsigned int nb, string &symbol,
1193  int &n1, float &v1, int &n2, float &v2, int &n3, float &v3)
1194 {
1195  string sbuf=string(buf);
1196  std::vector<std::string> v=SplitString(sbuf);
1197  if(nb==1)
1198  {//First atom "C 1"
1199  if(v.size()>0) symbol=v[0];
1200  else symbol=string(buf).substr(0,2);
1201  return;
1202  }
1203  if(nb==2)
1204  {//Second atom "C 1 1.450"
1205  if(v.size()==3)
1206  {
1207  symbol=v[0];
1208  n1=(unsigned int) atoi(v[1].c_str());
1209  v1=string2floatC(v[2]);
1210  }
1211  else
1212  {
1213  symbol=string(buf).substr(0,2);
1214  n1=(int) atoi(string(buf).substr(2,3).c_str());
1215  v1=string2floatC(string(buf).substr(5,6));
1216  }
1217  cout<<"ReadFHLine():"<<buf<<"#"<<symbol<<"/"<<n1<<"/"<<v1<<";"<<v.size()<<endl;
1218  VFN_DEBUG_MESSAGE("ReadFHLine():#"<<symbol<<"/"<<n1<<"/"<<v1,10);
1219  return;
1220  }
1221  if(nb==3)
1222  {//Third atom "C 2 1.450 1 119.995"
1223  if(v.size()==5)
1224  {
1225  symbol=v[0];
1226  n1=(int) atoi(v[1].c_str());
1227  v1=string2floatC(v[2]);
1228  n2=(int) atoi(v[3].c_str());
1229  v2=string2floatC(v[4]);
1230  }
1231  else
1232  {
1233  symbol=sbuf.substr(0,2);
1234  n1=(int) atoi(sbuf.substr(2,3).c_str());
1235  v1=string2floatC(sbuf.substr(5,6));
1236  n2=(int) atoi(sbuf.substr(11,3).c_str());
1237  v2=string2floatC(sbuf.substr(14,8));
1238  }
1239  VFN_DEBUG_MESSAGE("ReadFHLine():#"<<symbol<<"/"<<n1<<"/"<<v1<<"/"<<n2<<"/"<<v2,10);
1240  return;
1241  }
1242  //Remaining atoms
1243  if(v.size()==7)
1244  {
1245  symbol=v[0];
1246  n1=(int) atoi(v[1].c_str());
1247  v1=string2floatC(v[2]);
1248  n2=(int) atoi(v[3].c_str());
1249  v2=string2floatC(v[4]);
1250  n3=(int) atoi(v[5].c_str());
1251  v3=string2floatC(v[6]);
1252  }
1253  else
1254  {
1255  // sscanf ignores whitespace characters, so cannot be used ?? "H 601 1.100600 120.000599 180.0"
1256  //sscanf(buf,"%2s%3d%6f%3d%8f%3d%6f",symb,&n1,&v1,&n2,&v2,&n3,&v3);
1257  //symbol=string(symb);
1258  symbol=sbuf.substr(0,2);
1259  n1=(int) atoi(sbuf.substr(2,3).c_str());
1260  v1=string2floatC(sbuf.substr(5,6));
1261  n2=(int) atoi(sbuf.substr(11,3).c_str());
1262  v2=string2floatC(sbuf.substr(14,8));
1263  n3=(int) atoi(sbuf.substr(22,3).c_str());
1264  v3=string2floatC(sbuf.substr(25,6));
1265  }
1266  VFN_DEBUG_MESSAGE("ReadFHLine():#"<<symbol<<"/"<<n1<<"/"<<v1<<"/"<<n2<<"/"<<v2<<"/"<<n3<<"/"<<v3,10);
1267 }
1268 
1269 void ZScatterer::ImportFenskeHallZMatrix(istream &is,bool named)
1270 {
1271  char buf[101];
1272  // Get read of "KEYWORD GO HERE", just in case...
1273  {
1274  const char c=is.peek();
1275  if ( (c < '0') || (c > '9') )
1276  {
1277  cout<<"ZScatterer::ImportFenskeHallZMatrix()"
1278  <<":getting rid of first line..."<<endl;
1279  is.getline(buf,100);
1280  }
1281  }
1282  int nbAtoms=0;
1283  is >> nbAtoms;
1284  char c;
1285  while(!isgraph(is.peek())) is.get(c); // go to end of line
1286 
1287  // 17
1288  //C 1
1289  //N 1 1.465
1290  //C 2 1.366 1 119.987
1291  //N 3 1.321 2 120.030 1 6.0
1292  //C 4 1.355 3 119.982 2 6.8
1293  //N 5 1.136 4 180.000 3 46.3
1294  //N 3 1.366 2 120.022 1 186.0
1295  //C 7 1.466 3 119.988 2 354.9
1296  // ...
1297  string symbol, atomName,bondAtomName,angleAtomName,dihedAtomName,junk;
1298  int bondAtom=0,angleAtom=0,dihedAtom=0;
1299  float bond=0,angle=0,dihed=0;
1300  int scattPow;
1301  //first
1302  if(named)
1303  {
1304  is >> atomName >> symbol >> junk;
1305  VFN_DEBUG_MESSAGE("ZScatterer::ImportFenskeHallZMatrix():#"<<2<<",name:"<<atomName
1306  <<", bond: "<<bondAtomName<<", angle: "<<angleAtomName<<", dihed: "<<dihedAtomName,10);
1307  }
1308  else
1309  {
1310  is.getline(buf,100);
1311  ReadFHLine(buf,1,symbol,bondAtom,bond,angleAtom,angle,dihedAtom,dihed);
1312  }
1313  {
1314  scattPow=mpCryst->GetScatteringPowerRegistry().Find
1315  (symbol,"ScatteringPowerAtom",true);
1316  if(scattPow==-1)
1317  {
1318  cout<<"Scattering power"<<symbol<<"not found, creating it..."<<endl;
1319  mpCryst->AddScatteringPower(new ScatteringPowerAtom(symbol,symbol));
1320  }
1321  scattPow=mpCryst->GetScatteringPowerRegistry().Find
1322  (symbol,"ScatteringPowerAtom");
1323  if(named)
1324  this->AddAtom(atomName,
1325  &(mpCryst->GetScatteringPowerRegistry().GetObj(scattPow)),
1326  0,0,
1327  0,0,
1328  0,0);
1329  else
1330  {
1331  sprintf(buf,"%d",1);
1332  this->AddAtom(symbol+(string)buf,
1333  &(mpCryst->GetScatteringPowerRegistry().GetObj(scattPow)),
1334  0,0,
1335  0,0,
1336  0,0);
1337  }
1338  }
1339  //second
1340  if(named)
1341  {
1342  is >> atomName >> symbol >> bondAtomName >> bond;
1343  VFN_DEBUG_MESSAGE("ZScatterer::ImportFenskeHallZMatrix():#"<<2<<",name:"<<atomName
1344  <<", bond: "<<bondAtomName<<", angle: "<<angleAtomName<<", dihed: "<<dihedAtomName,10);
1345  }
1346  else
1347  {
1348  is.getline(buf,100);
1349  ReadFHLine(buf,2,symbol,bondAtom,bond,angleAtom,angle,dihedAtom,dihed);
1350  }
1351 
1352  {
1353  scattPow=mpCryst->GetScatteringPowerRegistry().Find
1354  (symbol,"ScatteringPowerAtom",true);
1355  if(scattPow==-1)
1356  {
1357  cout<<"Scattering power"<<symbol<<"not found, creating it..."<<endl;
1358  mpCryst->AddScatteringPower(new ScatteringPowerAtom(symbol,symbol));
1359  }
1360  scattPow=mpCryst->GetScatteringPowerRegistry().Find
1361  (symbol,"ScatteringPowerAtom");
1362  if(named)
1363  {
1364  bondAtom=mZAtomRegistry.Find(bondAtomName);
1365  if(bondAtom<0)
1366  throw ObjCrystException(string("ZScatterer::ImportFenskeHallZMatrix:")
1367  +string("when adding atom ")+atomName+string(": could not find atom: ")
1368  +bondAtomName);
1369  this->AddAtom(atomName,
1370  &(mpCryst->GetScatteringPowerRegistry().GetObj(scattPow)),
1371  bondAtom,bond,
1372  0,0,
1373  0,0);
1374  }
1375  else
1376  {
1377  sprintf(buf,"%d",2);
1378  this->AddAtom(symbol+(string)buf,
1379  &(mpCryst->GetScatteringPowerRegistry().GetObj(scattPow)),
1380  bondAtom-1,bond,
1381  0,0,
1382  0,0);
1383  }
1384  }
1385  //third
1386  if(named)
1387  {
1388  is >> atomName >> symbol >>
1389  bondAtomName >> bond >>
1390  angleAtomName >> angle;
1391  VFN_DEBUG_MESSAGE("ZScatterer::ImportFenskeHallZMatrix():#"<<2<<",name:"<<atomName
1392  <<", bond: "<<bondAtomName<<", angle: "<<angleAtomName<<", dihed: "<<dihedAtomName,10);
1393  }
1394  else
1395  {
1396  is.getline(buf,100);
1397  ReadFHLine(buf,3,symbol,bondAtom,bond,angleAtom,angle,dihedAtom,dihed);
1398  }
1399 
1400  {
1401  scattPow=mpCryst->GetScatteringPowerRegistry().Find
1402  (symbol,"ScatteringPowerAtom",true);
1403  if(scattPow==-1)
1404  {
1405  cout<<"Scattering power"<<symbol<<"not found, creating it..."<<endl;
1406  mpCryst->AddScatteringPower(new ScatteringPowerAtom(symbol,symbol));
1407  }
1408  scattPow=mpCryst->GetScatteringPowerRegistry().Find
1409  (symbol,"ScatteringPowerAtom");
1410  if(named)
1411  {
1412  VFN_DEBUG_MESSAGE("ZScatterer::ImportFenskeHallZMatrix():#"<<3<<",name:"<<atomName
1413  <<", bond: "<<bondAtomName<<", angle: "<<angleAtomName<<", dihed: "<<dihedAtomName,10);
1414  bondAtom=mZAtomRegistry.Find(bondAtomName);
1415  if(bondAtom<0)
1416  throw ObjCrystException(string("ZScatterer::ImportFenskeHallZMatrix:")
1417  +string("when adding atom ")+atomName+string(": could not find atom: ")
1418  +bondAtomName);
1419  angleAtom=mZAtomRegistry.Find(angleAtomName);
1420  if(angleAtom<0)
1421  throw ObjCrystException(string("ZScatterer::ImportFenskeHallZMatrix:")
1422  +string("when adding atom ")+atomName+string(": could not find atom: ")
1423  +angleAtomName);
1424  this->AddAtom(atomName,
1425  &(mpCryst->GetScatteringPowerRegistry().GetObj(scattPow)),
1426  bondAtom,bond,
1427  angleAtom,angle*DEG2RAD,
1428  0,0);
1429  }
1430  else
1431  {
1432  sprintf(buf,"%d",3);
1433  this->AddAtom(symbol+(string)buf,
1434  &(mpCryst->GetScatteringPowerRegistry().GetObj(scattPow)),
1435  bondAtom-1,bond,
1436  angleAtom-1,angle*DEG2RAD,
1437  0,0);
1438  }
1439  }
1440  for(int i=3;i<nbAtoms;i++)
1441  {
1442  if(named)
1443  {
1444  is >> atomName >> symbol >>
1445  bondAtomName >> bond >>
1446  angleAtomName >> angle >>
1447  dihedAtomName >> dihed;
1448  VFN_DEBUG_MESSAGE("ZScatterer::ImportFenskeHallZMatrix():#"<<i<<",name:"<<atomName
1449  <<", bond: "<<bondAtomName<<", angle: "<<angleAtomName<<", dihed: "<<dihedAtomName,10);
1450  }
1451  else
1452  {
1453  is.getline(buf,100);
1454  ReadFHLine(buf,i+1,symbol,bondAtom,bond,angleAtom,angle,dihedAtom,dihed);
1455  }
1456  {
1457  scattPow=mpCryst->GetScatteringPowerRegistry().Find
1458  (symbol,"ScatteringPowerAtom",true);
1459  if(scattPow==-1)
1460  {
1461  cout<<"Scattering power"<<symbol<<"not found, creating it..."<<endl;
1462  mpCryst->AddScatteringPower(new ScatteringPowerAtom(symbol,symbol));
1463  }
1464  scattPow=mpCryst->GetScatteringPowerRegistry().Find
1465  (symbol,"ScatteringPowerAtom");
1466  if(named)
1467  {
1468  bondAtom=mZAtomRegistry.Find(bondAtomName);
1469  angleAtom=mZAtomRegistry.Find(angleAtomName);
1470  dihedAtom=mZAtomRegistry.Find(dihedAtomName);
1471  VFN_DEBUG_MESSAGE("ZScatterer::ImportFenskeHallZMatrix():#"<<i<<",name:"<<atomName
1472  <<", bond: "<<bondAtomName<<", angle: "<<angleAtomName<<", dihed: "<<dihedAtomName,10);
1473  if(bondAtom<0)
1474  throw ObjCrystException(string("ZScatterer::ImportFenskeHallZMatrix:")
1475  +string("when adding atom ")+atomName+string(": could not find atom: ")
1476  +bondAtomName);
1477  if(bondAtom<0)
1478  throw ObjCrystException(string("ZScatterer::ImportFenskeHallZMatrix:")
1479  +string("when adding atom ")+atomName+string(": could not find atom: ")
1480  +angleAtomName);
1481  if(dihedAtom<0)
1482  throw ObjCrystException(string("ZScatterer::ImportFenskeHallZMatrix:")
1483  +string("when adding atom ")+atomName+string(": could not find atom: ")
1484  +dihedAtomName);
1485  this->AddAtom(atomName,
1486  &(mpCryst->GetScatteringPowerRegistry().GetObj(scattPow)),
1487  bondAtom,bond,
1488  angleAtom,angle*DEG2RAD,
1489  dihedAtom,dihed*DEG2RAD);
1490  }
1491  else
1492  {
1493  sprintf(buf,"%d",i+1);
1494  this->AddAtom(symbol+(string)buf,
1495  &(mpCryst->GetScatteringPowerRegistry().GetObj(scattPow)),
1496  bondAtom-1,bond,
1497  angleAtom-1,angle*DEG2RAD,
1498  dihedAtom-1,dihed*DEG2RAD);
1499  }
1500  }
1501  }
1502  this->SetLimitsRelative(gpRefParTypeScattConformBondLength,-.03,.03);
1503  this->SetLimitsRelative(gpRefParTypeScattConformBondAngle,-.01,.01);
1504  this->SetLimitsRelative(gpRefParTypeScattConformDihedAngle,-.01,.01);
1505 }
1506 
1507 void ZScatterer::ExportFenskeHallZMatrix(ostream &os)
1508 {
1509  if(mNbAtom<1) return;
1510  os << mNbAtom<<endl;
1511  os << this->GetZAtomRegistry().GetObj(0).mpScattPow->GetName()
1512  << " 1"<<endl;
1513  if(mNbAtom<2) return;
1514  os << this->GetZAtomRegistry().GetObj(1).mpScattPow->GetName()
1515  << " "<<this->GetZBondAtom(1)+1<< " "<<this->GetZBondLength(1)
1516  <<endl;
1517  if(mNbAtom<3) return;
1518  os << this->GetZAtomRegistry().GetObj(2).mpScattPow->GetName()
1519  << " "<<this->GetZBondAtom(2)+1 << " "<<this->GetZBondLength(2)
1520  << " "<<this->GetZAngleAtom(2)+1<< " "<<this->GetZAngle(2)*RAD2DEG
1521  <<endl;
1522  for(int i=3;i<mNbAtom;i++)
1523  os << this->GetZAtomRegistry().GetObj(i).mpScattPow->GetName()
1524  << " "<<this->GetZBondAtom(i)+1 << " "<<this->GetZBondLength(i)
1525  << " "<<this->GetZAngleAtom(i)+1<< " "<<this->GetZAngle(i)*RAD2DEG
1526  << " "<<this->GetZDihedralAngleAtom(i)+1<< " "<<this->GetZDihedralAngle(i)*RAD2DEG
1527  <<endl;
1528 }
1529 
1530 void ZScatterer::SetCenterAtomIndex(const unsigned int ix){mCenterAtomIndex=ix;}
1531 unsigned int ZScatterer::GetCenterAtomIndex()const{return mCenterAtomIndex;}
1532 
1533 std::size_t ZScatterer::size() const
1534 {
1535  return mZAtomRegistry.size();
1536 }
1537 
1538 vector<ZAtom*>::const_iterator ZScatterer::begin() const
1539 {
1540  return mZAtomRegistry.begin();
1541 }
1542 
1543 vector<ZAtom*>::const_iterator ZScatterer::end() const
1544 {
1545  return mZAtomRegistry.end();
1546 }
1547 
1548 void ZScatterer::GlobalOptRandomMove(const REAL mutationAmplitude,
1549  const RefParType *type)
1550 {
1551  if(mRandomMoveIsDone) return;
1552  VFN_DEBUG_ENTRY("ZScatterer::GlobalOptRandomMove()",3)
1553  TAU_PROFILE("ZScatterer::GlobalOptRandomMove()","void ()",TAU_DEFAULT);
1554  // give a 2% chance of either moving a single atom, or move
1555  // all atoms before a given torsion angle.
1556  // Only try this if there are more than 10 atoms (else it's not worth the speed cost)
1557  if((mNbAtom>=10) && ((rand()/(REAL)RAND_MAX)<.02)
1558  && (gpRefParTypeScattConform->IsDescendantFromOrSameAs(type)))//.01
1559  {
1560  TAU_PROFILE_TIMER(timer1,\
1561  "ZScatterer::GlobalOptRandomMoveSmart1(prepare ref par & mutate)"\
1562  ,"", TAU_FIELD);
1563  TAU_PROFILE_TIMER(timer2,\
1564  "ZScatterer::GlobalOptRandomMoveSmart2(optimize if necessary)"\
1565  ,"", TAU_FIELD);
1566  TAU_PROFILE_START(timer1);
1567  // Do we have *any* dihedral angle to really move ?
1568  CrystVector_long dihed(mNbAtom);
1569  dihed=0;
1570  int nbDihed=0;
1571  RefinablePar *par;
1572  dihed(nbDihed++)=2;//This is the Psi angle, in fact. Should Chi be added, too ?
1573  for(int i=3;i<mNbAtom;i++)
1574  {
1575  par=&(this->GetPar(&(mZAtomRegistry.GetObj(i).mDihed)));
1576  if( !(par->IsFixed()) ) //&& !(par->IsLimited())
1577  dihed(nbDihed++)=i;
1578  }
1579  if(nbDihed<2) //Can't play :-(
1580  {
1581  mRandomMoveIsDone=false;
1582  this->RefinableObj::GlobalOptRandomMove(mutationAmplitude);
1583  TAU_PROFILE_STOP(timer1);
1584  VFN_DEBUG_EXIT("ZScatterer::GlobalOptRandomMove():End",3)
1585  return;
1586  }
1587  // Build mpZMoveMinimizer object
1588  if(0==mpZMoveMinimizer)
1589  {
1590  mpZMoveMinimizer=new ZMoveMinimizer(*this);
1591  // copy all parameters (and not reference them, we will change the fix status..
1592  // we could remove all popu parameters...
1593  for(int i=0; i<this->GetNbPar();i++) mpZMoveMinimizer->AddPar(this->GetPar(i));
1594  }
1595  // Pick one to move and get the relevant parameter
1596  // (maybe we should random-move also the associated bond lengths an angles,
1597  // but for now we'll concentrate on dihedral (torsion) angles.
1598  const int atom=dihed((int) (rand()/((REAL)RAND_MAX+1)*nbDihed));
1599  //cout<<endl;
1600  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove(): Changing atom #"<<atom ,3)
1601  if(atom==2)
1602  par=&(this->GetPar(&mPsi));
1603  else
1604  par=&(this->GetPar(&(mZAtomRegistry.GetObj(atom).mDihed)));
1605  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove(): initial value:"<<par->GetHumanValue() ,3)
1606  // Record the current conformation
1607  mpZMoveMinimizer->RecordConformation();
1608  // Set up
1609  const int moveType= rand()%3;
1610  mpZMoveMinimizer->FixAllPar();
1611  REAL x0,y0,z0;
1612  //cout << " Move Type:"<<moveType<<endl;
1613  CrystVector_REAL weight(mNbAtom);
1614  switch(moveType)
1615  {// :TODO: rather build a connectivity table, to include more atoms
1616  case 0:// (0) Try to move only one atom
1617  {
1618  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove():Move only one atom",3)
1619  weight=1;
1620  weight(atom)=0;
1621  for(int i=0;i<mNbAtom;i++)
1622  if( (i==this->GetZBondAtom(atom))
1623  ||(i==this->GetZAngleAtom(atom))
1624  ||(i==this->GetZDihedralAngleAtom(atom))
1625  ||(atom==this->GetZBondAtom(i))
1626  ||(atom==this->GetZAngleAtom(i))
1627  ||(atom==this->GetZDihedralAngleAtom(i)))
1628  {
1629  if(!(this->GetPar(&(mZAtomRegistry.GetObj(i).mBondLength)).IsFixed()))
1630  mpZMoveMinimizer->GetPar(&(mZAtomRegistry.GetObj(i).mBondLength)).SetIsFixed(false);
1631  if(!(this->GetPar(&(mZAtomRegistry.GetObj(i).mAngle)).IsFixed()))
1632  mpZMoveMinimizer->GetPar(&(mZAtomRegistry.GetObj(i).mAngle)).SetIsFixed(false);
1633  if(!(this->GetPar(&(mZAtomRegistry.GetObj(i).mDihed)).IsFixed()))
1634  mpZMoveMinimizer->GetPar(&(mZAtomRegistry.GetObj(i).mDihed)).SetIsFixed(false);
1635  }
1636  if(!(this->GetPar(&mPhi).IsFixed()))
1637  mpZMoveMinimizer->GetPar(&mPhi).SetIsFixed(false);
1638  if(!(this->GetPar(&mChi).IsFixed()))
1639  mpZMoveMinimizer->GetPar(&mChi).SetIsFixed(false);
1640  if( !(this->GetPar(&mPsi).IsFixed()) && (atom!=2))
1641  mpZMoveMinimizer->GetPar(&mPsi).SetIsFixed(false);
1642  if(!(this->GetPar(&mXYZ(0)).IsFixed()))
1643  mpZMoveMinimizer->GetPar(&mXYZ(0)).SetIsFixed(false);
1644  if(!(this->GetPar(&mXYZ(1)).IsFixed()))
1645  mpZMoveMinimizer->GetPar(&mXYZ(1)).SetIsFixed(false);
1646  if(!(this->GetPar(&mXYZ(2)).IsFixed()))
1647  mpZMoveMinimizer->GetPar(&mXYZ(2)).SetIsFixed(false);
1648  break;
1649  }
1650  case 1:// (1) Try to move the atoms *before* the rotated bond
1651  {
1652  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove():Move before the rotated bond",3)
1653  const int atom1=this->GetZBondAtom(atom);
1654  const int atom2=this->GetZAngleAtom(atom);
1655  weight=0;
1656  weight(atom1)=1;
1657  for(int i=0;i<mNbAtom;i++) if(weight(this->GetZBondAtom(i))>.1) weight(i)=1;
1658  weight(atom2)=1;
1659  for(int i=0;i<mNbAtom;i++)
1660  if( (atom2 ==this->GetZAngleAtom(i))
1661  &&(atom1 !=this->GetZBondAtom(i))
1662  &&(i != atom) )
1663  {
1664  if(!(this->GetPar(&(mZAtomRegistry.GetObj(i).mBondLength)).IsFixed()))
1665  mpZMoveMinimizer->GetPar(&(mZAtomRegistry.GetObj(i).mBondLength)).SetIsFixed(false);
1666  if(!(this->GetPar(&(mZAtomRegistry.GetObj(i).mAngle)).IsFixed()))
1667  mpZMoveMinimizer->GetPar(&(mZAtomRegistry.GetObj(i).mAngle)).SetIsFixed(false);
1668  if(!(this->GetPar(&(mZAtomRegistry.GetObj(i).mDihed)).IsFixed()))
1669  mpZMoveMinimizer->GetPar(&(mZAtomRegistry.GetObj(i).mDihed)).SetIsFixed(false);
1670  }
1671  if(!(this->GetPar(&mPhi).IsFixed()))
1672  mpZMoveMinimizer->GetPar(&mPhi).SetIsFixed(false);
1673  if(!(this->GetPar(&mChi).IsFixed()))
1674  mpZMoveMinimizer->GetPar(&mChi).SetIsFixed(false);
1675  if( !(this->GetPar(&mPsi).IsFixed()) && (atom!=2))
1676  mpZMoveMinimizer->GetPar(&mPsi).SetIsFixed(false);
1677 
1678  if(!(this->GetPar(&mXYZ(0)).IsFixed()))
1679  mpZMoveMinimizer->GetPar(&mXYZ(0)).SetIsFixed(false);
1680  if(!(this->GetPar(&mXYZ(1)).IsFixed()))
1681  mpZMoveMinimizer->GetPar(&mXYZ(1)).SetIsFixed(false);
1682  if(!(this->GetPar(&mXYZ(2)).IsFixed()))
1683  mpZMoveMinimizer->GetPar(&mXYZ(2)).SetIsFixed(false);
1684 
1685  if(!(this->GetPar(&(mZAtomRegistry.GetObj(atom).mBondLength)).IsFixed()))
1686  mpZMoveMinimizer->GetPar(&(mZAtomRegistry.GetObj(atom).mBondLength)).SetIsFixed(false);
1687  if(!(this->GetPar(&(mZAtomRegistry.GetObj(atom).mAngle)).IsFixed()))
1688  mpZMoveMinimizer->GetPar(&(mZAtomRegistry.GetObj(atom).mAngle)).SetIsFixed(false);
1689 
1690  break;
1691  }
1692  case 2:// (2) Try to move the atoms *after* the rotated bond
1693  {
1694  if(mCenterAtomIndex>=atom)
1695  {
1696  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove():Move after the bond (translate)",3)
1697  x0=this->GetXCoord()(0);
1698  y0=this->GetYCoord()(0);
1699  z0=this->GetZCoord()(0);
1700  }
1701  else
1702  {
1703  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove():Move after the bond(nothing to do)",3)
1704  }
1705  break;
1706  }
1707  }
1708  // Move it, and with some probability use flipping to some
1709  // not-so-random angles., and then minimize the conformation change
1710  mpZMoveMinimizer->SetZAtomWeight(weight);
1711  REAL change;
1712  if( (rand()%5)==0)
1713  {
1714  switch(rand()%5)
1715  {
1716  case 0: change=-120*DEG2RAD;break;
1717  case 1: change= -90*DEG2RAD;break;
1718  case 2: change= 90*DEG2RAD;break;
1719  case 3: change= 120*DEG2RAD;break;
1720  default:change= 180*DEG2RAD;break;
1721  }
1722  }
1723  else
1724  {
1725  change= par->GetGlobalOptimStep()
1726  *2*(rand()/(REAL)RAND_MAX-0.5)*mutationAmplitude*16;
1727  }
1728  TAU_PROFILE_STOP(timer1);
1729  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove(): mutation:"<<change*RAD2DEG,3)
1730  par->Mutate(change);
1731  if(2==moveType)
1732  {
1733  if(mCenterAtomIndex>=atom)
1734  {
1735  this->UpdateCoordinates();
1736  x0 -= this->GetXCoord()(0);
1737  y0 -= this->GetYCoord()(0);
1738  z0 -= this->GetZCoord()(0);
1739  mpCryst->OrthonormalToFractionalCoords(x0,y0,z0);
1740  this->GetPar(mXYZ.data()).Mutate(x0);
1741  this->GetPar(mXYZ.data()+1).Mutate(y0);
1742  this->GetPar(mXYZ.data()+2).Mutate(z0);
1743  }
1744  }
1745  else
1746  {
1747  const REAL tmp=mpZMoveMinimizer->GetLogLikelihood();
1748  if(tmp>.05)
1749  {
1750  TAU_PROFILE_START(timer2);
1751  if(tmp<1) mpZMoveMinimizer->MinimizeChange(100);
1752  else if(tmp<5) mpZMoveMinimizer->MinimizeChange(200);
1753  else mpZMoveMinimizer->MinimizeChange(500);
1754  TAU_PROFILE_STOP(timer2);
1755  }
1756  }
1757  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove(): final value:"<<par->GetHumanValue(),3)
1758  //
1759  }
1760  #if 0
1761  {
1762  // find unfixed, dihedral angles to play with //unlimited?
1763  CrystVector_long dihed(mNbAtom);
1764  dihed=0;
1765  int nbDihed=0;
1766  RefinablePar *par;
1767  dihed(nbDihed++)=2;
1768  for(int i=3;i<mNbAtom;i++)
1769  {
1770  par=&(this->GetPar(&(mZAtomRegistry.GetObj(i).mDihed)));
1771  if( !(par->IsFixed()) ) //&& !(par->IsLimited())
1772  dihed(nbDihed++)=i;
1773  }
1774  if(nbDihed<2) //Can't play :-(
1775  this->RefinableObj::GlobalOptRandomMove(mutationAmplitude);
1776  // Pick one
1777  const int atom=dihed((int) (rand()/((REAL)RAND_MAX+1)*nbDihed));
1778  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove(): "<<FormatHorizVector<long>(dihed) ,10)
1779  VFN_DEBUG_MESSAGE("ZScatterer::GlobalOptRandomMove(): Changing atom #"<<atom ,10)
1780  if(atom==2)
1781  par=&(this->GetPar(&mPsi));
1782  else
1783  par=&(this->GetPar(&(mZAtomRegistry.GetObj(atom).mDihed)));
1784  // Get the old value
1785  const REAL old=par->GetValue();
1786  // Move it, with a max amplitude 8x greater than usual
1787  if( (rand()/(REAL)RAND_MAX)<.1)
1788  {// give some probability to use certain angles: -120,-90,90,120,180
1789  switch(rand()%5)
1790  {
1791  case 0: par->Mutate(-120*!DEG2RAD);break;
1792  case 1: par->Mutate( -90*!DEG2RAD);break;
1793  case 2: par->Mutate( 90*!DEG2RAD);break;
1794  case 3: par->Mutate( 120*!DEG2RAD);break;
1795  default:par->Mutate( 180*!DEG2RAD);break;
1796  }
1797  }
1798  else
1799  par->Mutate( par->GetGlobalOptimStep()
1800  *2*(rand()/(REAL)RAND_MAX-0.5)*mutationAmplitude*8);
1801  const REAL change=mZAtomRegistry.GetObj(atom).GetZDihedralAngle()-old;
1802  // Now move all atoms using this changed bond as a reference
1803  //const int atom2= mZAtomRegistry.GetObj(atom).GetZAngleAtom();
1804  for(int i=atom;i<mNbAtom;i++)
1805  //if( (mZAtomRegistry.GetObj(i).GetZBondAtom()==atom)
1806  // &&(mZAtomRegistry.GetObj(i).GetZAngleAtom()==atom2))
1807  // this->GetPar(&(mZAtomRegistry.GetObj(i).mDihed)).Mutate(-change);
1808  if(mZAtomRegistry.GetObj(i).GetZAngleAtom()==atom)
1809  this->GetPar(&(mZAtomRegistry.GetObj(i).mDihed)).Mutate(-change);
1810  //cout <<"ZScatterer::GlobalOptRandomMove:"<<nbDihed
1811  // <<" "<<atom
1812  // <<" "<<atom2
1813  // <<" "<<rand()
1814  // <<endl
1815  // <<" "<<FormatHorizVector<long>(dihed,4)
1816  // <<endl;
1817  }
1818  #endif
1819  else
1820  {
1821  mRandomMoveIsDone=false;
1822  this->RefinableObj::GlobalOptRandomMove(mutationAmplitude,type);
1823  }
1824  mRandomMoveIsDone=true;
1825  VFN_DEBUG_EXIT("ZScatterer::GlobalOptRandomMove():End",3)
1826 }
1827 
1828 void ZScatterer::UpdateCoordinates() const
1829 {
1830  if(mClockCoord>mClockScatterer) return;
1831 
1832  VFN_DEBUG_ENTRY("ZScatterer::UpdateCoordinates():"<<this->GetName(),3)
1833  TAU_PROFILE("ZScatterer::UpdateCoordinates()","void ()",TAU_DEFAULT);
1834  //if(0==mNbAtom) throw ObjCrystException("ZScatterer::Update() No Atoms in Scatterer !");
1835  if(0==mNbAtom) return;
1836 
1837  {
1838  CrystMatrix_REAL phiMatrix(3,3),chiMatrix(3,3),psiMatrix(3,3);
1839  phiMatrix= cos(mPhi) , -sin(mPhi) , 0,
1840  sin(mPhi) , cos(mPhi) , 0,
1841  0 ,0 ,1;
1842 
1843  chiMatrix= cos(mChi) ,0 ,-sin(mChi),
1844  0 ,1 ,0,
1845  sin(mChi) ,0 ,cos(mChi);
1846 
1847  psiMatrix= 1 , 0 , 0,
1848  0 ,cos(mPsi) ,-sin(mPsi),
1849  0 ,sin(mPsi) ,cos(mPsi);
1850 
1851  mPhiChiPsiMatrix=product(chiMatrix,product(phiMatrix,psiMatrix));
1852  //cout << phiMatrix <<endl<< chiMatrix <<endl<< psiMatrix <<endl<<mPhiChiPsiMatrix<<endl;
1853  }
1854 
1855  mXCoord.resize(mNbAtom);
1856  mYCoord.resize(mNbAtom);
1857  mZCoord.resize(mNbAtom);
1858 
1859  // Atom 0
1860  mXCoord(0)=0.;
1861  mYCoord(0)=0.;
1862  mZCoord(0)=0.;
1863  VFN_DEBUG_MESSAGE("->Atom #0:"<<mXCoord(0)<<" : "<<mYCoord(0)<<" : "<<mZCoord(0),1)
1864 
1865  if(mNbAtom>1)
1866  {// Atom 1
1867  mXCoord(1)=GetZBondLength(1);
1868  mYCoord(1)=0.;
1869  mZCoord(1)=0.;
1870  VFN_DEBUG_MESSAGE("->Atom #1:"<<mXCoord(1)<<" : "<<mYCoord(1)<<" : "<<mZCoord(1),1)
1871  }
1872  if(mNbAtom>2)
1873  {// Atom 2
1874  if(0==GetZBondAtom(2)) //Linked with Atom 1
1875  mXCoord(2)=GetZBondLength(2)*cos(GetZAngle(2));
1876  else //Linked with Atom 1
1877  mXCoord(2)=mXCoord(1)-GetZBondLength(2)*cos(GetZAngle(2));
1878  mYCoord(2)=GetZBondLength(2)*sin(GetZAngle(2));
1879  mZCoord(2)=0.;
1880  VFN_DEBUG_MESSAGE("->Atom #2:"<<mXCoord(2)<<" : "<<mYCoord(2)<<" : "<<mZCoord(2),1)
1881  }
1882  for(int i=1;i<3;i++)// Global rotation of scatterer
1883  {
1884  if(mNbAtom==i)break;
1885  const REAL x=mXCoord(i);
1886  const REAL y=mYCoord(i);
1887  const REAL z=mZCoord(i);
1888  mXCoord(i)=mPhiChiPsiMatrix(0,0)*x+mPhiChiPsiMatrix(0,1)*y+mPhiChiPsiMatrix(0,2)*z;
1889  mYCoord(i)=mPhiChiPsiMatrix(1,0)*x+mPhiChiPsiMatrix(1,1)*y+mPhiChiPsiMatrix(1,2)*z;
1890  mZCoord(i)=mPhiChiPsiMatrix(2,0)*x+mPhiChiPsiMatrix(2,1)*y+mPhiChiPsiMatrix(2,2)*z;
1891  }
1892  if(mNbAtom>3)
1893  {
1894  REAL xa,ya,za,xb,yb,zb,xd,yd,zd,cosph,sinph,costh,sinth,coskh,sinkh,cosa,sina;
1895  REAL xpd,ypd,zpd,xqd,yqd,zqd;
1896  REAL rbc,xyb,yza,tmp,xpa,ypa,zqa;
1897  int na,nb,nc;
1898  bool flag;
1899  REAL dist,angle,dihed;
1900  for(int i=3;i<mNbAtom;i++)
1901  {
1902  na=GetZBondAtom(i);
1903  nb=GetZAngleAtom(i);
1904  nc=GetZDihedralAngleAtom(i);
1905  dist=GetZBondLength(i);
1906  angle=GetZAngle(i);
1907  dihed=GetZDihedralAngle(i);
1908 
1909  xb = mXCoord(nb) - mXCoord(na);
1910  yb = mYCoord(nb) - mYCoord(na);
1911  zb = mZCoord(nb) - mZCoord(na);
1912 
1913  rbc= sqrt(xb*xb + yb*yb + zb*zb);
1914  if(rbc<1e-5)
1915  {
1916  throw ObjCrystException("ZScatterer::UpdateCoordinates(): two atoms ("
1917  +mZAtomRegistry.GetObj(na).GetName()
1918  +" and "+ mZAtomRegistry.GetObj(nb).GetName()
1919  +") have the same coordinates (d<1e-5): aborting.");
1920  }
1921  rbc=1./rbc;
1922 
1923  cosa = cos(angle);
1924  sina = sin(angle);
1925 
1926  if( fabs(cosa) >= 0.999999 )
1927  { // Colinear
1928  mXCoord(i)=mXCoord(na)+cosa*dist*rbc*xb;
1929  mYCoord(i)=mYCoord(na)+cosa*dist*rbc*yb;
1930  mZCoord(i)=mZCoord(na)+cosa*dist*rbc*zb;
1931  VFN_DEBUG_MESSAGE("->Atom #"<<i<<":"<<mXCoord(i)<<" : "<<mYCoord(i)<<" : " <<mZCoord(i)<<"(colinear)",1)
1932  }
1933  else
1934  {
1935  xa = mXCoord(nc) - mXCoord(na);
1936  ya = mYCoord(nc) - mYCoord(na);
1937  za = mZCoord(nc) - mZCoord(na);
1938  xd = dist*cosa;
1939  yd = dist*sina*cos(dihed);
1940  zd = -dist*sina*sin(dihed);
1941 
1942  xyb = sqrt(xb*xb + yb*yb);
1943  if( xyb < 0.1 )
1944  { // Rotate about y-axis
1945  tmp = za; za = -xa; xa = tmp;
1946  tmp = zb; zb = -xb; xb = tmp;
1947  xyb = sqrt(xb*xb + yb*yb);
1948  flag = true;
1949  }
1950  else flag = false;
1951 
1952  costh = xb/xyb;
1953  sinth = yb/xyb;
1954  xpa = costh*xa + sinth*ya;
1955  ypa = costh*ya - sinth*xa;
1956  sinph = zb*rbc;
1957  cosph = sqrt(1.0 - sinph*sinph);
1958  zqa = cosph*za - sinph*xpa;
1959  yza = sqrt(ypa*ypa + zqa*zqa);
1960 
1961  //cout <<endl<<ypa<<" "<<zqa<<" "<<yza<<" "<<endl;
1962  if( yza > 1e-5 )
1963  {
1964  coskh = ypa/yza;
1965  sinkh = zqa/yza;
1966  ypd = coskh*yd - sinkh*zd;
1967  zpd = coskh*zd + sinkh*yd;
1968  }
1969  else
1970  {
1971  ypd = yd;
1972  zpd = zd;
1973  }
1974 
1975  xpd = cosph*xd - sinph*zpd;
1976  zqd = cosph*zpd + sinph*xd;
1977  xqd = costh*xpd - sinth*ypd;
1978  yqd = costh*ypd + sinth*xpd;
1979 
1980  if( true==flag )
1981  { // Rotate about y-axis ?
1982  mXCoord(i)=mXCoord(na) - zqd;
1983  mYCoord(i)=mYCoord(na) + yqd;
1984  mZCoord(i)=mZCoord(na) + xqd;
1985  } else
1986  {
1987  mXCoord(i)=mXCoord(na) + xqd;
1988  mYCoord(i)=mYCoord(na) + yqd;
1989  mZCoord(i)=mZCoord(na) + zqd;
1990  }
1991  VFN_DEBUG_MESSAGE("->Atom #"<<i<<":"<<mXCoord(i)<<" : "<<mYCoord(i)<<" : " <<mZCoord(i),1)
1992  }
1993  }
1994  }
1995  //shift atom around Central atom
1996  REAL x,y,z;
1997  x=this->GetX();
1998  y=this->GetY();
1999  z=this->GetZ();
2000  mpCryst->FractionalToOrthonormalCoords(x,y,z);
2001  const REAL x0=x-mXCoord(mCenterAtomIndex);
2002  const REAL y0=y-mYCoord(mCenterAtomIndex);
2003  const REAL z0=z-mZCoord(mCenterAtomIndex);
2004  for(int i=0;i<mNbAtom;i++)
2005  {
2006  mXCoord(i) += x0;
2007  mYCoord(i) += y0;
2008  mZCoord(i) += z0;
2009  }
2010  mClockCoord.Click();
2011  VFN_DEBUG_EXIT("ZScatterer::UpdateCoordinates()"<<this->GetName(),3)
2012 }
2013 
2014 void ZScatterer::UpdateScattCompList() const
2015 {
2016  VFN_DEBUG_ENTRY("ZScatterer::UpdateScattCompList()"<<this->GetName(),3)
2017  this->UpdateCoordinates();
2018  if( (mClockScattCompList>mClockCoord)
2019  &&(mClockScattCompList>mpCryst->GetClockLatticePar())) return;
2020 
2021  if(true==mUseGlobalScattPow)
2022  {
2023  mScattCompList(0).mX=mXYZ(0);
2024  mScattCompList(0).mY=mXYZ(1);
2025  mScattCompList(0).mZ=mXYZ(2);
2026  mScattCompList(0).mOccupancy=mOccupancy;
2027  mScattCompList(0).mpScattPow=mpGlobalScattPow;
2028 
2029  // Update central atom for Display.
2030  //mXCoord(mCenterAtomIndex)=this->GetX();
2031  //mYCoord(mCenterAtomIndex)=this->GetY();
2032  //mZCoord(mCenterAtomIndex)=this->GetZ();
2033  //mpCryst->FractionalToOrthonormalCoords(mXCoord(mCenterAtomIndex),
2034  // mXCoord(mCenterAtomIndex),
2035  // mXCoord(mCenterAtomIndex));
2036 
2037  mClockScattCompList.Click();
2038  VFN_DEBUG_MESSAGE("ZScatterer::UpdateScattCompList()->Global Scatterer:End",3)
2039  return;
2040  }
2041 
2042  long j=0;
2043  REAL x,y,z;
2044  VFN_DEBUG_MESSAGE("ZScatterer::UpdateScattCompList(bool):Finishing"<<mNbAtom<<","<<mNbDummyAtom,3)
2045  for(long i=0;i<mNbAtom;i++)
2046  {
2047  if(0!=mZAtomRegistry.GetObj(i).GetScatteringPower())
2048  {
2049  x=mXCoord(i);
2050  y=mYCoord(i);
2051  z=mZCoord(i);
2052  mpCryst->OrthonormalToFractionalCoords(x,y,z);
2053  mScattCompList(j ).mX=x;
2054  mScattCompList(j ).mY=y;
2055  mScattCompList(j ).mZ=z;
2056  mScattCompList(j ).mpScattPow=mZAtomRegistry.GetObj(i).GetScatteringPower();
2057  mScattCompList(j++).mOccupancy=mZAtomRegistry.GetObj(i).GetOccupancy()*mOccupancy;
2058  }
2059  }
2060  #ifdef __DEBUG__
2061  if(gVFNDebugMessageLevel<3) mScattCompList.Print();
2062  #endif
2063  mClockScattCompList.Click();
2064  VFN_DEBUG_EXIT("ZScatterer::UpdateScattCompList()"<<this->GetName(),3)
2065 }
2066 
2067 void ZScatterer::InitRefParList()
2068 {
2069  VFN_DEBUG_MESSAGE("ZScatterer::InitRefParList():"<<this->GetName(),5)
2070  //throw ObjCrystException("ZScatterer::InitRefParList() not implemented! "+this->GetName());
2071  //:TODO:
2072  this->ResetParList();
2073  {
2074  RefinablePar tmp("x",&mXYZ(0),0.,1.,
2075  gpRefParTypeScattTranslX,
2076  REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,1.,1.);
2077  tmp.AssignClock(mClockScatterer);
2078  this->AddPar(tmp);
2079  }
2080  {
2081  RefinablePar tmp("y",&mXYZ(1),0,1,
2082  gpRefParTypeScattTranslY,
2083  REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,1.,1.);
2084  tmp.AssignClock(mClockScatterer);
2085  this->AddPar(tmp);
2086  }
2087  {
2088  RefinablePar tmp("z",&mXYZ(2),0,1,
2089  gpRefParTypeScattTranslZ,
2090  REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,1.,1.);
2091  tmp.AssignClock(mClockScatterer);
2092  this->AddPar(tmp);
2093  }
2094  {
2095  RefinablePar tmp("Occupancy",&mOccupancy,0,1,
2096  gpRefParTypeScattOccup,
2097  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,1.,1.);
2098  tmp.AssignClock(mClockScatterer);
2099  this->AddPar(tmp);
2100  }
2101  if(false==mUseGlobalScattPow)
2102  {
2103  {
2104  RefinablePar tmp("Phi",&mPhi,0,2*M_PI,
2105  gpRefParTypeScattOrient,
2106  REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,RAD2DEG,2*M_PI);
2107  tmp.AssignClock(mClockScatterer);
2108  this->AddPar(tmp);
2109  }
2110  {
2111  RefinablePar tmp("Chi",&mChi,0,2*M_PI,
2112  gpRefParTypeScattOrient,
2113  REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,RAD2DEG,2*M_PI);
2114  tmp.AssignClock(mClockScatterer);
2115  this->AddPar(tmp);
2116  }
2117  {
2118  RefinablePar tmp("Psi",&mPsi,0,2*M_PI,
2119  gpRefParTypeScattOrient,
2120  REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,RAD2DEG,2*M_PI);
2121  tmp.AssignClock(mClockScatterer);
2122  this->AddPar(tmp);
2123  }
2124  char buf [20];
2125  for(long i=0;i<mNbAtom;i++)
2126  {
2127  bool usedBond=true,usedAngle=true,usedDihed=true;
2128  if(i<1) usedBond=false;
2129  if(i<2) usedAngle=false;
2130  if(i<3) usedDihed=false;
2131  {
2132  sprintf(buf,"%d-%d",(int)i,(int)(mZAtomRegistry.GetObj(i).GetZBondAtom()));
2133  RefinablePar tmp((string)"Length"+(string)buf,
2134  &(mZAtomRegistry.GetObj(i).mBondLength),
2135  mZAtomRegistry.GetObj(i).mBondLength*.9,
2136  mZAtomRegistry.GetObj(i).mBondLength*1.1,
2137  gpRefParTypeScattConformBondLength,
2138  REFPAR_DERIV_STEP_ABSOLUTE,true,false,usedBond,false,1.);
2139  tmp.AssignClock(mClockScatterer);
2140  this->AddPar(tmp);
2141  }
2142  {
2143  sprintf(buf,"%d-%d-%d",(int)i,(int)(mZAtomRegistry.GetObj(i).GetZBondAtom()),
2144  (int)(mZAtomRegistry.GetObj(i).GetZAngleAtom()));
2145  RefinablePar tmp("Angle"+(string)buf,
2146  &(mZAtomRegistry.GetObj(i).mAngle),0,2*M_PI,
2147  gpRefParTypeScattConformBondAngle,
2148  REFPAR_DERIV_STEP_ABSOLUTE,false,false,usedAngle,true,RAD2DEG,2*M_PI);
2149  tmp.AssignClock(mClockScatterer);
2150  this->AddPar(tmp);
2151  }
2152  {
2153  sprintf(buf,"%d-%d-%d-%d",(int)i,(int)(mZAtomRegistry.GetObj(i).GetZBondAtom()),
2154  (int)(mZAtomRegistry.GetObj(i).GetZAngleAtom()),
2155  (int)(mZAtomRegistry.GetObj(i).GetZDihedralAngleAtom()));
2156  RefinablePar tmp("Dihed"+(string)buf,
2157  &(mZAtomRegistry.GetObj(i).mDihed),0,2*M_PI,
2158  gpRefParTypeScattConformDihedAngle,
2159  REFPAR_DERIV_STEP_ABSOLUTE,false,false,usedDihed,true,RAD2DEG,2*M_PI);
2160  tmp.AssignClock(mClockScatterer);
2161  this->AddPar(tmp);
2162  }
2163  if(0!=mZAtomRegistry.GetObj(i).GetScatteringPower())
2164  {//fixed by default
2165  sprintf(buf,"%d",(int)i);
2166  RefinablePar tmp("Occupancy"+(string)buf,
2167  &(mZAtomRegistry.GetObj(i).mOccupancy),0,1,
2168  gpRefParTypeScattOccup,
2169  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,1.,1.);
2170  tmp.AssignClock(mClockScatterer);
2171  this->AddPar(tmp);
2172  }
2173  }
2174  }//if(false==mUseGlobalScatteringPower)
2175 }
2176 
2177 #ifdef __WX__CRYST__
2178 WXCrystObjBasic* ZScatterer::WXCreate(wxWindow* parent)
2179 {
2180  //:TODO: Check mpWXCrystObj==0
2181  mpWXCrystObj=new WXZScatterer(parent,this);
2182  return mpWXCrystObj;
2183 }
2184 #endif
2185 
2186 //######################################################################
2187 //
2188 // ZPolyhedron
2189 //
2190 //
2191 //######################################################################
2192 ZPolyhedron::ZPolyhedron( const RegularPolyhedraType type,Crystal &cryst,
2193  const REAL x, const REAL y, const REAL z,
2194  const string &name, const ScatteringPower *centralAtomSymbol,
2195  const ScatteringPower *periphAtomSymbol,const REAL centralPeriphDist,
2196  const REAL ligandPopu,
2197  const REAL phi, const REAL chi, const REAL psi):
2198 ZScatterer(name,cryst,x,y,z,phi,chi,psi),mPolyhedraType(type)
2199 {
2200  VFN_DEBUG_MESSAGE("ZPolyhedron::ZPolyhedron(..)",5)
2201  // Additioning string and char arrays takes a huge lot of mem (gcc 2.95.2)
2202  const string name_=name+"_";
2203  const string name_central=name_+centralAtomSymbol->GetName();
2204  const string name_periph=name_+periphAtomSymbol->GetName();
2205  switch(mPolyhedraType)
2206  {
2207  case TETRAHEDRON :
2208  {
2209  REAL ang=2*asin(sqrt(2./3.));
2210  this ->AddAtom (name_central, centralAtomSymbol,
2211  0,0.,
2212  0,0.,
2213  0,0.,
2214  1.);
2215  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2216  0,centralPeriphDist,
2217  0,0.,
2218  0,0.,
2219  1.);
2220  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2221  0,centralPeriphDist,
2222  1,ang,
2223  0,0.,
2224  1.);
2225  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2226  0,centralPeriphDist,
2227  1,ang,
2228  2, M_PI*2./3.,
2229  1.);
2230  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2231  0,centralPeriphDist,
2232  1,ang,
2233  2,-M_PI*2./3.,
2234  1.);
2235  m3DDisplayIndex.resize(4,3);
2236  m3DDisplayIndex= 1,2,3,
2237  1,2,4,
2238  1,3,4,
2239  2,3,4;
2240  break;
2241  }
2242  case OCTAHEDRON:
2243  {
2244  this ->AddAtom (name_central, centralAtomSymbol,
2245  0,0.,
2246  0,0.,
2247  0,0.,
2248  1.);
2249  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2250  0,centralPeriphDist,
2251  0,0.,
2252  0,0.,
2253  1.);
2254  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2255  0,centralPeriphDist,
2256  1,M_PI/2,
2257  0,0.,
2258  1.);
2259  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2260  0,centralPeriphDist,
2261  1,M_PI/2,
2262  2,M_PI/2,
2263  1.);
2264  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2265  0,centralPeriphDist,
2266  1,M_PI/2,
2267  2,-M_PI/2,
2268  1.);
2269  this ->AddAtom (name_periph+"5",periphAtomSymbol,
2270  0,centralPeriphDist,
2271  1,M_PI/2,
2272  2,M_PI,
2273  1.);
2274  this ->AddAtom (name_periph+"6",periphAtomSymbol,
2275  0,centralPeriphDist,
2276  1,M_PI,
2277  2,0.,
2278  1.);
2279  m3DDisplayIndex.resize(8,3);
2280  m3DDisplayIndex= 1,2,3,
2281  1,2,4,
2282  1,5,3,
2283  1,5,4,
2284  6,2,3,
2285  6,2,4,
2286  6,5,3,
2287  6,5,4;
2288  break;
2289  }
2290  case SQUARE_PLANE:
2291  {
2292  this ->AddAtom (name_central, centralAtomSymbol,
2293  0,0.,
2294  0,0.,
2295  0,0.,
2296  1.);
2297  this ->AddAtom (name+"_X",0, //Dummy atom on top
2298  0,1.,
2299  0,0.,
2300  0,0.,
2301  1.);
2302  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2303  0,centralPeriphDist,
2304  1,M_PI/2,
2305  0,0.,
2306  1.);
2307  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2308  0,centralPeriphDist,
2309  1,M_PI/2,
2310  2,M_PI/2,
2311  1.);
2312  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2313  0,centralPeriphDist,
2314  1,M_PI/2,
2315  2,M_PI,
2316  1.);
2317  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2318  0,centralPeriphDist,
2319  1,M_PI/2,
2320  2,-M_PI/2,
2321  1.);
2322  //:TODO: GL Display with squares
2323  //m3DDisplayIndex.resize(2,3);
2324  //m3DDisplayIndex= 2,4,2,
2325  // 2,4,5;
2326  break;
2327  }
2328  case CUBE:
2329  {
2330  this ->AddAtom (name_central, centralAtomSymbol,
2331  0,0.,
2332  0,0.,
2333  0,0.,
2334  1.);
2335  this ->AddAtom (name+"_X",0, //Dummy atom in middle of face
2336  0,1.,
2337  0,0.,
2338  0,0.,
2339  1.);
2340  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2341  0,centralPeriphDist,
2342  1,M_PI/4.,
2343  0,0.,
2344  1.);
2345  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2346  0,centralPeriphDist,
2347  1,M_PI/4.,
2348  2,M_PI/2.,
2349  1.);
2350  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2351  0,centralPeriphDist,
2352  1,M_PI/4.,
2353  2,M_PI,
2354  1.);
2355  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2356  0,centralPeriphDist,
2357  1,M_PI/4.,
2358  2,-M_PI/2.,
2359  1.);
2360  this ->AddAtom (name_periph+"5",periphAtomSymbol,
2361  0,centralPeriphDist,
2362  1,M_PI*3./4.,
2363  0,0.,
2364  1.);
2365  this ->AddAtom (name_periph+"6",periphAtomSymbol,
2366  0,centralPeriphDist,
2367  1,M_PI*3./4.,
2368  2,M_PI/2.,
2369  1.);
2370  this ->AddAtom (name_periph+"7",periphAtomSymbol,
2371  0,centralPeriphDist,
2372  1,M_PI*3./4.,
2373  2,M_PI,
2374  1.);
2375  this ->AddAtom (name_periph+"8",periphAtomSymbol,
2376  0,centralPeriphDist,
2377  1,M_PI*3./4.,
2378  2,M_PI*3./2.,
2379  1.);
2380  break;
2381  }
2382  case ANTIPRISM_TETRAGONAL:
2383  {
2384  this ->AddAtom (name_central, centralAtomSymbol,
2385  0,0.,
2386  0,0.,
2387  0,0.,
2388  1.);
2389  this ->AddAtom (name+"_X",0, //Dummy atom in middle of face
2390  0,1.,
2391  0,0.,
2392  0,0.,
2393  1.);
2394  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2395  0,centralPeriphDist,
2396  1,M_PI/4.,
2397  0,0.,
2398  1.);
2399  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2400  0,centralPeriphDist,
2401  1,M_PI/4.,
2402  2,M_PI/2.,
2403  1.);
2404  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2405  0,centralPeriphDist,
2406  1,M_PI/4.,
2407  2,M_PI,
2408  1.);
2409  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2410  0,centralPeriphDist,
2411  1,M_PI/4.,
2412  2,-M_PI/2.,
2413  1.);
2414  this ->AddAtom (name_periph+"5",periphAtomSymbol,
2415  0,centralPeriphDist,
2416  1,M_PI*3./4.,
2417  0,M_PI/4.,
2418  1.);
2419  this ->AddAtom (name_periph+"6",periphAtomSymbol,
2420  0,centralPeriphDist,
2421  1,M_PI*3./4.,
2422  2,M_PI*3./4.,
2423  1.);
2424  this ->AddAtom (name_periph+"7",periphAtomSymbol,
2425  0,centralPeriphDist,
2426  1,M_PI*3./4.,
2427  2,M_PI*5./4.,
2428  1.);
2429  this ->AddAtom (name_periph+"8",periphAtomSymbol,
2430  0,centralPeriphDist,
2431  1,M_PI*3./4.,
2432  2,M_PI*7./4.,
2433  1.);
2434  break;
2435  }
2436  case PRISM_TETRAGONAL_MONOCAP:
2437  {
2438  this ->AddAtom (name_central, centralAtomSymbol,
2439  0,0.,
2440  0,0.,
2441  0,0.,
2442  1.);
2443  this ->AddAtom (name_periph+"0",periphAtomSymbol,
2444  0,centralPeriphDist,
2445  0,0.,
2446  0,0.,
2447  1.);
2448  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2449  0,centralPeriphDist,
2450  1,70*DEG2RAD,
2451  0,0.,
2452  1.);
2453  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2454  0,centralPeriphDist,
2455  1,70*DEG2RAD,
2456  2,M_PI/2.,
2457  1.);
2458  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2459  0,centralPeriphDist,
2460  1,70*DEG2RAD,
2461  2,M_PI,
2462  1.);
2463  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2464  0,centralPeriphDist,
2465  1,70*DEG2RAD,
2466  2,-M_PI/2.,
2467  1.);
2468  this ->AddAtom (name_periph+"5",periphAtomSymbol,
2469  0,centralPeriphDist,
2470  1,145*DEG2RAD,
2471  0,0.,
2472  1.);
2473  this ->AddAtom (name_periph+"6",periphAtomSymbol,
2474  0,centralPeriphDist,
2475  1,145*DEG2RAD,
2476  2,M_PI/2.,
2477  1.);
2478  this ->AddAtom (name_periph+"7",periphAtomSymbol,
2479  0,centralPeriphDist,
2480  1,145*DEG2RAD,
2481  2,M_PI,
2482  1.);
2483  this ->AddAtom (name_periph+"8",periphAtomSymbol,
2484  0,centralPeriphDist,
2485  1,145*DEG2RAD,
2486  2,M_PI*3./2.,
2487  1.);
2488  break;
2489  }
2490  case PRISM_TETRAGONAL_DICAP:
2491  {
2492  this ->AddAtom (name_central, centralAtomSymbol,
2493  0,0.,
2494  0,0.,
2495  0,0.,
2496  1.);
2497  this ->AddAtom (name_periph+"0",periphAtomSymbol,
2498  0,centralPeriphDist,
2499  0,0.,
2500  0,0.,
2501  1.);
2502  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2503  0,centralPeriphDist,
2504  1,60*DEG2RAD,
2505  0,0.,
2506  1.);
2507  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2508  0,centralPeriphDist,
2509  1,60*DEG2RAD,
2510  2,M_PI/2.,
2511  1.);
2512  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2513  0,centralPeriphDist,
2514  1,60*DEG2RAD,
2515  2,M_PI,
2516  1.);
2517  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2518  0,centralPeriphDist,
2519  1,60*DEG2RAD,
2520  2,-M_PI/2.,
2521  1.);
2522  this ->AddAtom (name_periph+"5",periphAtomSymbol,
2523  0,centralPeriphDist,
2524  1,120*DEG2RAD,
2525  0,0.,
2526  1.);
2527  this ->AddAtom (name_periph+"6",periphAtomSymbol,
2528  0,centralPeriphDist,
2529  1,120*DEG2RAD,
2530  2,M_PI/2.,
2531  1.);
2532  this ->AddAtom (name_periph+"7",periphAtomSymbol,
2533  0,centralPeriphDist,
2534  1,120*DEG2RAD,
2535  2,M_PI,
2536  1.);
2537  this ->AddAtom (name_periph+"8",periphAtomSymbol,
2538  0,centralPeriphDist,
2539  1,120*DEG2RAD,
2540  2,M_PI*3./2.,
2541  1.);
2542  this ->AddAtom (name_periph+"8",periphAtomSymbol,
2543  0,centralPeriphDist,
2544  1,M_PI,
2545  2,0.,
2546  1.);
2547  break;
2548  }
2549  case PRISM_TRIGONAL:
2550  {
2551  const REAL ang=55.*DEG2RAD;
2552  const REAL ang2=120.*DEG2RAD;
2553  this ->AddAtom (name_central, centralAtomSymbol,
2554  0,0.,
2555  0,0.,
2556  0,0.,
2557  1.);
2558  this ->AddAtom (name+"_X",0, //Dummy atom in middle of top face
2559  0,1.,
2560  0,0.,
2561  0,0.,
2562  0.);
2563  this ->AddAtom (name_periph+"0",periphAtomSymbol,
2564  0,centralPeriphDist,
2565  1,ang,
2566  0,0.,
2567  1.);
2568  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2569  0,centralPeriphDist,
2570  1,ang,
2571  2,ang2,
2572  1.);
2573  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2574  0,centralPeriphDist,
2575  1,ang,
2576  2,-ang2,
2577  1.);
2578  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2579  0,centralPeriphDist,
2580  1,M_PI-ang,
2581  0,0.,
2582  1.);
2583  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2584  0,centralPeriphDist,
2585  1,M_PI-ang,
2586  2,ang2,
2587  1.);
2588  this ->AddAtom (name_periph+"5",periphAtomSymbol,
2589  0,centralPeriphDist,
2590  1,M_PI-ang,
2591  2,-ang2,
2592  1.);
2593  break;
2594  }
2595  case PRISM_TRIGONAL_TRICAPPED:
2596  {
2597  const REAL ang=55.*DEG2RAD;
2598  const REAL ang2=120.*DEG2RAD;
2599  this ->AddAtom (name_central, centralAtomSymbol,
2600  0,0.,
2601  0,0.,
2602  0,0.,
2603  1.);
2604  this ->AddAtom (name+"_X",0, //Dummy atom in middle of top face
2605  0,1.,
2606  0,0.,
2607  0,0.,
2608  0.);
2609  this ->AddAtom (name_periph+"0",periphAtomSymbol,
2610  0,centralPeriphDist,
2611  1,ang,
2612  0,0.,
2613  1.);
2614  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2615  0,centralPeriphDist,
2616  1,ang,
2617  2,ang2,
2618  1.);
2619  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2620  0,centralPeriphDist,
2621  1,ang,
2622  2,-ang2,
2623  1.);
2624  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2625  0,centralPeriphDist,
2626  1,M_PI-ang,
2627  0,0.,
2628  1.);
2629  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2630  0,centralPeriphDist,
2631  1,M_PI-ang,
2632  2,ang2,
2633  1.);
2634  this ->AddAtom (name_periph+"5",periphAtomSymbol,
2635  0,centralPeriphDist,
2636  1,M_PI-ang,
2637  2,-ang2,
2638  1.);
2639  this ->AddAtom (name_periph+"6",periphAtomSymbol,
2640  0,centralPeriphDist,
2641  1,M_PI/2.,
2642  2,M_PI,
2643  1.);
2644  this ->AddAtom (name_periph+"7",periphAtomSymbol,
2645  0,centralPeriphDist,
2646  1,M_PI/2.,
2647  2,M_PI/3.,
2648  1.);
2649  this ->AddAtom (name_periph+"8",periphAtomSymbol,
2650  0,centralPeriphDist,
2651  1,M_PI/2.,
2652  2,-M_PI/3.,
2653  1.);
2654  break;
2655  }
2656  case ICOSAHEDRON:
2657  {
2658  const REAL ang=acos(sqrt(.2));
2659  const REAL ang2=M_PI*2./5.;
2660  this ->AddAtom (name_central, centralAtomSymbol,
2661  0,0.,
2662  0,0.,
2663  0,0.,
2664  1.);
2665  this ->AddAtom (name_periph+"0",periphAtomSymbol,
2666  0,centralPeriphDist,
2667  0,0.,
2668  0,0.,
2669  1.);
2670  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2671  0,centralPeriphDist,
2672  1,ang,
2673  1,0.,
2674  1.);
2675  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2676  0,centralPeriphDist,
2677  1,ang,
2678  2,ang2,
2679  1.);
2680  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2681  0,centralPeriphDist,
2682  1,ang,
2683  2,2*ang2,
2684  1.);
2685  this ->AddAtom (name_periph+"4",periphAtomSymbol,
2686  0,centralPeriphDist,
2687  1,ang,
2688  2,-2*ang2,
2689  1.);
2690  this ->AddAtom (name_periph+"5",periphAtomSymbol,
2691  0,centralPeriphDist,
2692  1,ang,
2693  2,-ang2,
2694  1.);
2695  this ->AddAtom (name_periph+"6",periphAtomSymbol,
2696  0,centralPeriphDist,
2697  1,M_PI,
2698  0,0.,
2699  1.);
2700  this ->AddAtom (name_periph+"7",periphAtomSymbol,
2701  0,centralPeriphDist,
2702  7,ang,
2703  3,M_PI,
2704  1.);
2705  this ->AddAtom (name_periph+"8",periphAtomSymbol,
2706  0,centralPeriphDist,
2707  7,ang,
2708  8,ang2,
2709  1.);
2710  this ->AddAtom (name_periph+"9",periphAtomSymbol,
2711  0,centralPeriphDist,
2712  7,ang,
2713  8,2*ang2,
2714  1.);
2715  this ->AddAtom (name_periph+"10",periphAtomSymbol,
2716  0,centralPeriphDist,
2717  7,ang,
2718  8,-2*ang2,
2719  1.);
2720  this ->AddAtom (name_periph+"11",periphAtomSymbol,
2721  0,centralPeriphDist,
2722  7,ang,
2723  8,-ang2,
2724  1.);
2725  break;
2726  }
2727  case TRIANGLE_PLANE:
2728  {
2729  this ->AddAtom (name_central, centralAtomSymbol,
2730  0,0.,
2731  0,0.,
2732  0,0.,
2733  1.);
2734  this ->AddAtom (name+"_X",0, //Dummy atom on top
2735  0,1.,
2736  0,0.,
2737  0,0.,
2738  1.);
2739  this ->AddAtom (name_periph+"1",periphAtomSymbol,
2740  0,centralPeriphDist,
2741  1,M_PI/2,
2742  0,0.,
2743  1.);
2744  this ->AddAtom (name_periph+"2",periphAtomSymbol,
2745  0,centralPeriphDist,
2746  1,M_PI/2,
2747  2,2*M_PI/3,
2748  1.);
2749  this ->AddAtom (name_periph+"3",periphAtomSymbol,
2750  0,centralPeriphDist,
2751  1,M_PI/2,
2752  2,-2*M_PI/3,
2753  1.);
2754  //:TODO: GL Display with squares
2755  //m3DDisplayIndex.resize(2,3);
2756  //m3DDisplayIndex= 2,4,2,
2757  // 2,4,5;
2758  break;
2759  }
2760  default : throw ObjCrystException("ZPolyhedron::ZPolyhedron():Unknown Polyhedra type !");
2761  }
2762  //We want to keep an approximate geometry
2763  //this->RefinableObj::Print();
2764  this->SetLimitsRelative(gpRefParTypeScattConformBondLength,-.2,.2);
2765  this->SetLimitsRelative(gpRefParTypeScattConformBondAngle ,-.2,.2);
2766  this->SetLimitsRelative(gpRefParTypeScattConformDihedAngle,-.2,.2);
2767  //this->RefinableObj::Print();
2768 
2769  VFN_DEBUG_MESSAGE("ZPolyhedron::ZPolyhedron():End:"<<mName<<")",5)
2770 }
2771 
2772 ZPolyhedron::ZPolyhedron(const ZPolyhedron &old):
2773 ZScatterer(old){}
2774 
2775 ZPolyhedron* ZPolyhedron::CreateCopy() const
2776 {
2777  VFN_DEBUG_MESSAGE("ZPolyhedron::CreateCopy()"<<mName<<")",5)
2778  return new ZPolyhedron(*this);
2779 }
2780 
2781 //######################################################################
2782 //
2783 // GLOBAL SCATTERING POWER
2784 //
2785 //######################################################################
2786 
2787 GlobalScatteringPower::GlobalScatteringPower():mpZScatterer(0)
2788 {
2789  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GlobalScatteringPower():"<<mName,5)
2790 }
2791 
2792 GlobalScatteringPower::GlobalScatteringPower(const ZScatterer &scatt):mpZScatterer(0)
2793 {
2794  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GlobalScatteringPower(&scatt):"<<mName,5)
2795  this->Init(scatt);
2796 }
2797 
2798 GlobalScatteringPower::GlobalScatteringPower(const GlobalScatteringPower& old):mpZScatterer(0)
2799 {
2800  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GlobalScatteringPower(&old):"<<mName,5)
2801  this->Init(*(old.mpZScatterer));
2802 }
2803 
2804 GlobalScatteringPower::~GlobalScatteringPower()
2805 {
2806  VFN_DEBUG_MESSAGE("GlobalScatteringPower::~GlobalScatteringPower():"<<mName,5)
2807  if(mpZScatterer!=0) delete mpZScatterer;
2808 }
2809 
2810 // Disable the base-class function.
2811 void GlobalScatteringPower::Init()
2812 {
2813  // This should be never called.
2814  abort();
2815 }
2816 
2817 void GlobalScatteringPower::Init(const ZScatterer &scatt)
2818 {
2819  this->SetName(scatt.GetName()+(string)"_Global");
2820  VFN_DEBUG_MESSAGE("GlobalScatteringPower::Init(&Scatt):"<<mName,5)
2821  if(mpZScatterer!=0) delete mpZScatterer;
2822  mpZScatterer=new ZScatterer(scatt);
2823  mpZScatterer->SetUseGlobalScatteringPower(false);
2824  mpZScatterer->SetName(scatt.GetName()+"_GlobalCopy");
2825 
2826  //Set the DynPopCorrIndex to the sum of the DynPopCorrIndexs (eg the sum of atomic numbers)
2827  mDynPopCorrIndex=0;
2828 
2829  const ScatteringComponentList* tmp=&(mpZScatterer->GetScatteringComponentList());
2830  for(long i=0;i<tmp->GetNbComponent();i++)
2831  mDynPopCorrIndex += (*tmp)(i).mpScattPow->GetDynPopCorrIndex();
2832 }
2833 
2834 CrystVector_REAL GlobalScatteringPower::
2835  GetScatteringFactor(const ScatteringData &data,
2836  const int spgSymPosIndex) const
2837 {
2838  // Here comes the hard work
2839  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GetScatteringFactor():"<<mName,10)
2840  TAU_PROFILE("GlobalScatteringPower::GetScatteringFactor()","void ()",TAU_DEFAULT);
2841 
2842  // copy both the scatterer and the DiffractionData object to determine
2843  // the average isotropic scattering power for each reflection
2844  ScatteringData* pData=data.CreateCopy();
2845  CrystVector_REAL sf(data.GetNbRefl()),rsf,isf;
2846  sf=0;
2847 
2848  Crystal cryst(data.GetCrystal().GetLatticePar(0),data.GetCrystal().GetLatticePar(1),
2849  data.GetCrystal().GetLatticePar(2),data.GetCrystal().GetLatticePar(3),
2850  data.GetCrystal().GetLatticePar(4),data.GetCrystal().GetLatticePar(5),"P1");
2851  cryst.AddScatterer(mpZScatterer);
2852  cryst.SetUseDynPopCorr(false);
2853  pData->SetCrystal(cryst);
2854  pData->SetName("GlobalScatteringPowerData!");
2855  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GetScatteringFactor():No DEBUG Messages",5)
2856  VFN_DEBUG_LOCAL_LEVEL(10)
2857 
2858  const long nbStep=4;//number of steps over 90 degrees for phi,chi psi
2859  REAL norm=0;
2860  for(int i=-nbStep;i<=nbStep;i++)
2861  {
2862  mpZScatterer->SetChi(i*M_PI/2/nbStep);
2863  for(int j=-nbStep;j<=nbStep;j++)
2864  {
2865  mpZScatterer->SetPhi(j*M_PI/2/nbStep);
2866  for(int k=-nbStep;k<=nbStep;k++)
2867  {
2868  //cout <<i<<","<<j<<","<<k<<endl;
2869  mpZScatterer->SetPsi(k*M_PI/2/nbStep);
2870 
2871  rsf=pData->GetFhklCalcReal();
2872  rsf*=rsf;
2873  isf=pData->GetFhklCalcImag();
2874  isf*=isf;
2875  rsf+=isf;
2876  rsf=sqrt(rsf);
2877 
2878  //correct for the solid angle (dChi*dPhi) corresponding to this orientation
2879  rsf*= cos(mpZScatterer->GetPhi());//:TODO: needs checking !!!
2880  norm += cos(mpZScatterer->GetPhi());
2881  sf += rsf;
2882  }
2883  }
2884  //cout << FormatHorizVector<REAL>(sf) <<endl<<norm<<endl;
2885  }
2886 
2887  VFN_DEBUG_LOCAL_LEVEL(-1)
2888  sf /= norm;
2889  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GetScatteringFactor():"<<mName<<":End.",10)
2890  delete pData;
2891  return sf;
2892 }
2893 REAL GlobalScatteringPower::GetForwardScatteringFactor(const RadiationType type) const
2894 {
2895  REAL sf=0;
2896  const ScatteringComponentList *pList=&(mpZScatterer->GetScatteringComponentList());
2897  for(int i=0;i<pList->GetNbComponent();i++)
2898  sf += (*pList)(i).mpScattPow->GetForwardScatteringFactor(type);
2899  return sf;
2900 }
2901 CrystVector_REAL GlobalScatteringPower::
2902  GetTemperatureFactor(const ScatteringData &data,
2903  const int spgSymPosIndex) const
2904 {
2905  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GetTemperatureFactor(data):"<<mName,5)
2906  CrystVector_REAL temp(data.GetNbRefl());
2907  temp=1.;
2908  return temp;
2909 }
2910 CrystMatrix_REAL GlobalScatteringPower::
2911  GetResonantScattFactReal(const ScatteringData &data,
2912  const int spgSymPosIndex) const
2913 {
2914  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GetResonantScattFactReal(data):"<<mName,5)
2915  CrystMatrix_REAL res(1,1);
2916  res=0.;
2917  return res;
2918 }
2919 CrystMatrix_REAL GlobalScatteringPower::
2920  GetResonantScattFactImag(const ScatteringData &data,
2921  const int spgSymPosIndex) const
2922 {
2923  VFN_DEBUG_MESSAGE("GlobalScatteringPower::GetResonantScattFactImag(data):"<<mName,5)
2924  CrystMatrix_REAL res(1,1);
2925  res=0.;
2926  return res;
2927 }
2928 REAL GlobalScatteringPower::GetRadius()const
2929 {
2930  //:TODO:
2931  return 3.;
2932 }
2933 void GlobalScatteringPower::InitRefParList()
2934 {
2935  //nothing to do, nothing to refine 8-))
2936 }
2937 
2938 }//namespace
ObjCryst
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++.
Definition: doc-main.h:25
ObjCryst::ReadFHLine
void ReadFHLine(const char *buf, const unsigned int nb, string &symbol, int &n1, float &v1, int &n2, float &v2, int &n3, float &v3)
Function to parse one line from a Fenske-Hall zmatrix file.
Definition: ZScatterer.cpp:1192
ObjCryst::RadiationType
RadiationType
Type of radiation used.
Definition: General.h:96
ObjCryst::string2floatC
float string2floatC(const string &s)
Function to convert a substring to a floating point value, imposing a C locale (using '.
Definition: ObjCryst/IO.cpp:59
ObjCryst::Crystal
Crystal class: Unit cell, spacegroup, scatterers.
Definition: Crystal.h:98
ObjCryst::Crystal::SetUseDynPopCorr
void SetUseDynPopCorr(const int use)
Set the use of dynamical population correction (Crystal::mUseDynPopCorr).
Definition: Crystal.cpp:859
ObjCryst::Crystal::AddScatteringPower
void AddScatteringPower(ScatteringPower *scattPow)
Add a ScatteringPower for this Crystal.
Definition: Crystal.cpp:241
ObjCryst::Crystal::GetScatteringPowerRegistry
ObjRegistry< ScatteringPower > & GetScatteringPowerRegistry()
Get the registry of ScatteringPower included in this Crystal.
Definition: Crystal.cpp:236
ObjCryst::Crystal::AddScatterer
void AddScatterer(Scatterer *scatt)
Add a scatterer to the crystal.
Definition: Crystal.cpp:188
ObjCryst::ObjCrystException
Exception class for ObjCryst++ library.
Definition: General.h:122
ObjCryst::CrystalPOVRayOptions
Class to store POV-Ray output options.
Definition: General.h:178
ObjCryst::Scatterer
Generic type of scatterer: can be an atom, or a more complex assembly of atoms.
Definition: Scatterer.h:131
ObjCryst::Scatterer::mpCryst
Crystal * mpCryst
The crystal in which the Scatterer is This is needed so that we can know which scattering powers are ...
Definition: Scatterer.h:304
ObjCryst::Scatterer::GetCrystal
const Crystal & GetCrystal() const
In which crystal is this Scatterer included ?
Definition: Scatterer.cpp:137
ObjCryst::Scatterer::mClockScattCompList
RefinableObjClock mClockScattCompList
Definition: Scatterer.h:297
ObjCryst::Scatterer::mXYZ
CrystVector_REAL mXYZ
coordinates of the scatterer (or of its center..)
Definition: Scatterer.h:283
ObjCryst::Scatterer::GetX
REAL GetX() const
X coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
Definition: Scatterer.cpp:103
ObjCryst::Scatterer::GetZ
REAL GetZ() const
Z coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
Definition: Scatterer.cpp:105
ObjCryst::Scatterer::mClockScatterer
RefinableObjClock mClockScatterer
Last time anything (number of atoms, positions, scattering power) was changed.
Definition: Scatterer.h:295
ObjCryst::Scatterer::SetCrystal
void SetCrystal(Crystal &)
Set the crystal in which is included this Scatterer.
Definition: Scatterer.cpp:136
ObjCryst::Scatterer::mOccupancy
REAL mOccupancy
Occupancy : 0 <= occ <= 1 For a multi-atom scatterer (polyhedron,..), this is the overall occupancy o...
Definition: Scatterer.h:289
ObjCryst::Scatterer::GetY
REAL GetY() const
Y coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
Definition: Scatterer.cpp:104
ObjCryst::ScatteringData
Class to compute structure factors for a set of reflections and a Crystal.
Definition: ScatteringData.h:320
ObjCryst::ScatteringData::SetCrystal
virtual void SetCrystal(Crystal &crystal)
Set the crystal for this experiment.
Definition: ScatteringData.cpp:654
ObjCryst::ScatteringData::CreateCopy
virtual ScatteringData * CreateCopy() const =0
So-called virtual copy constructor.
ObjCryst::ScatteringData::GetFhklCalcImag
const CrystVector_REAL & GetFhklCalcImag() const
Access to imaginary part of F(hkl)calc.
Definition: ScatteringData.cpp:791
ObjCryst::ScatteringData::GetNbRefl
long GetNbRefl() const
Return the number of reflections in this experiment.
Definition: ScatteringData.cpp:673
ObjCryst::ScatteringData::GetCrystal
const Crystal & GetCrystal() const
Const access to the data's crystal.
Definition: ScatteringData.cpp:667
ObjCryst::ScatteringData::GetFhklCalcReal
const CrystVector_REAL & GetFhklCalcReal() const
Access to real part of F(hkl)calc.
Definition: ScatteringData.cpp:783
ObjCryst::ScatteringPower
Abstract Base Class to describe the scattering power of any Scatterer component in a crystal.
Definition: ScatteringPower.h:111
ObjCryst::ScatteringPowerAtom
The Scattering Power for an Atom.
Definition: ScatteringPower.h:347
ObjCryst::ScatteringComponentList
list of scattering positions in a crystal, associated with the corresponding occupancy and a pointer ...
Definition: ScatteringPower.h:506
ObjCryst::ScatteringComponentList::GetNbComponent
long GetNbComponent() const
Number of components.
Definition: ScatteringPower.cpp:995
ObjCryst::ScatteringComponentList::Print
void Print() const
Print the list of Scattering components. For debugging.
Definition: ScatteringPower.cpp:1035
ObjCryst::UnitCell::GetLatticePar
CrystVector_REAL GetLatticePar() const
Lattice parameters (a,b,c,alpha,beta,gamma) as a 6-element vector in Angstroems and radians.
Definition: UnitCell.cpp:92
ObjCryst::UnitCell::OrthonormalToFractionalCoords
void OrthonormalToFractionalCoords(REAL &x, REAL &y, REAL &z) const
Get fractional cartesian coordinates for a set of (x,y,z) orthonormal coordinates.
Definition: UnitCell.cpp:273
ObjCryst::UnitCell::GetClockLatticePar
const RefinableObjClock & GetClockLatticePar() const
last time the Lattice parameters were changed
Definition: UnitCell.cpp:333
ObjCryst::UnitCell::FractionalToOrthonormalCoords
void FractionalToOrthonormalCoords(REAL &x, REAL &y, REAL &z) const
Get orthonormal cartesian coordinates for a set of (x,y,z) fractional coordinates.
Definition: UnitCell.cpp:263
ObjCryst::GlobalScatteringPower
Global Scattering Power.
Definition: ZScatterer.h:53
ObjCryst::GlobalScatteringPower::GetTemperatureFactor
virtual CrystVector_REAL GetTemperatureFactor(const ScatteringData &data, const int spgSymPosIndex=0) const
Get the temperature factor for all reflections of a given ScatteringData object.
Definition: ZScatterer.cpp:2902
ObjCryst::GlobalScatteringPower::GetForwardScatteringFactor
virtual REAL GetForwardScatteringFactor(const RadiationType) const
Get the scattering factor at (0,0,0).
Definition: ZScatterer.cpp:2893
ObjCryst::GlobalScatteringPower::GetResonantScattFactReal
virtual CrystMatrix_REAL GetResonantScattFactReal(const ScatteringData &data, const int spgSymPosIndex=0) const
Get the real part of the resonant scattering factor.
Definition: ZScatterer.cpp:2911
ObjCryst::GlobalScatteringPower::Init
void Init()
Initialization of the object, used by all constructors, and operator=.
Definition: ZScatterer.cpp:2811
ObjCryst::GlobalScatteringPower::GetResonantScattFactImag
virtual CrystMatrix_REAL GetResonantScattFactImag(const ScatteringData &data, const int spgSymPosIndex=0) const
Get the imaginary part of the resonant scattering factor.
Definition: ZScatterer.cpp:2920
ObjCryst::GlobalScatteringPower::GetRadius
virtual REAL GetRadius() const
Return the physical radius of this type of scatterer (for 3D display purposes).
Definition: ZScatterer.cpp:2928
ObjCryst::GlobalScatteringPower::GetScatteringFactor
virtual CrystVector_REAL GetScatteringFactor(const ScatteringData &data, const int spgSymPosIndex=0) const
Get the Scattering factor for all reflections of a given ScatteringData object.
Definition: ZScatterer.cpp:2835
ObjCryst::ZAtom
Class for individual atoms in a ZScatterer Object.
Definition: ZScatterer.h:87
ObjCryst::ZAtom::mBondLength
REAL mBondLength
Bond length, angle and dihedral angle.
Definition: ZScatterer.h:143
ObjCryst::ZMoveMinimizer
Class to minimize conformation changes for random moves.
Definition: ZScatterer.h:169
ObjCryst::ZMoveMinimizer::GetLogLikelihood
virtual REAL GetLogLikelihood() const
Get -log(likelihood) of the current configuration for the object.
Definition: ZScatterer.cpp:153
ObjCryst::ZScatterer
ZScatterer: the basic type of complex scatterers, where atom positions are defined using a standard "...
Definition: ZScatterer.h:191
ObjCryst::ZScatterer::ImportFenskeHallZMatrix
void ImportFenskeHallZMatrix(istream &is, bool named=false)
Import "Fenske-Hall" ZMatrix file (fhz in the babel program http://www.eyesopen.com/babel....
Definition: ZScatterer.cpp:1269
ObjCryst::ZScatterer::GetZAtomX
REAL GetZAtomX(const int i) const
Get the X fractionnal coordinate of atom i.
Definition: ZScatterer.cpp:438
ObjCryst::ZScatterer::SetZAngle
void SetZAngle(const int i, const REAL)
Access to the angle parameter, for the i-th row in the Z-Matrix.
Definition: ZScatterer.cpp:460
ObjCryst::ZScatterer::GlobalOptRandomMove
virtual void GlobalOptRandomMove(const REAL mutationAmplitude, const RefParType *type=gpRefParTypeObjCryst)
Make a random move of the current configuration.
Definition: ZScatterer.cpp:1548
ObjCryst::ZScatterer::SetCenterAtomIndex
void SetCenterAtomIndex(const unsigned int)
Set the index of the central atom (around which the rotation is made)
Definition: ZScatterer.cpp:1530
ObjCryst::ZScatterer::ZScatterer
ZScatterer(const string &name, Crystal &cryst, const REAL x=0., const REAL y=0., const REAL z=0., const REAL phi=0., const REAL chi=0., const REAL psi=0.)
ZScatterer constructor.
Definition: ZScatterer.cpp:214
ObjCryst::ZScatterer::size
std::size_t size() const
STL access to registry of Zatom size.
Definition: ZScatterer.cpp:1533
ObjCryst::ZScatterer::CreateCopy
virtual ZScatterer * CreateCopy() const
Definition: ZScatterer.cpp:306
ObjCryst::ZScatterer::GetZAtomRegistry
const ObjRegistry< ZAtom > & GetZAtomRegistry() const
Access to the registry of ZAtoms.
Definition: ZScatterer.cpp:466
ObjCryst::ZScatterer::mpGlobalScattPow
GlobalScatteringPower * mpGlobalScattPow
the global scattering power used, if mUseGlobalScattPow=true
Definition: ZScatterer.h:455
ObjCryst::ZScatterer::mPhi
REAL mPhi
Angles giving the orientation of the ZScatterer (stored in radian)
Definition: ZScatterer.h:438
ObjCryst::ZScatterer::GetZBondAtom
long GetZBondAtom(const int i) const
Index of the 1st atom used to define the i-th atom in the Z-Matrix (the one from which the bondlength...
Definition: ZScatterer.cpp:441
ObjCryst::ZScatterer::GetZAtomZ
REAL GetZAtomZ(const int i) const
Get the Z fractionnal coordinate of atom i.
Definition: ZScatterer.cpp:440
ObjCryst::ZScatterer::InitRefParList
virtual void InitRefParList()
Prepare refinable parameters for the scatterer object.
Definition: ZScatterer.cpp:2067
ObjCryst::ZScatterer::ExportFenskeHallZMatrix
void ExportFenskeHallZMatrix(ostream &os)
Export to Fenske-Hall ZMatrix file.
Definition: ZScatterer.cpp:1507
ObjCryst::ZScatterer::mClockCoord
RefinableObjClock mClockCoord
Last time the cartesian coordinates were computed.
Definition: ZScatterer.h:461
ObjCryst::ZScatterer::mCenterAtomIndex
long mCenterAtomIndex
Index of the atom used as a pivot (the scatterer is rotated around this atom).
Definition: ZScatterer.h:444
ObjCryst::ZScatterer::GetNbComponent
virtual int GetNbComponent() const
Number of components in the scatterer (eg number of point scatterers)
Definition: ZScatterer.cpp:405
ObjCryst::ZScatterer::mXCoord
CrystVector_REAL mXCoord
Storage for Cartesian coordinates.
Definition: ZScatterer.h:459
ObjCryst::ZScatterer::GetXCoord
const CrystVector_REAL & GetXCoord() const
Get the list of all ZAtom cartesian x coordinates.
Definition: ZScatterer.cpp:1130
ObjCryst::ZScatterer::GetPhi
REAL GetPhi() const
Access to phi parameter (overall orientation of the scatterer)
Definition: ZScatterer.cpp:430
ObjCryst::ZScatterer::SetZDihedralAngle
void SetZDihedralAngle(const int i, const REAL)
Access to the dihedral angle parameter, for the i-th row in the Z-Matrix.
Definition: ZScatterer.cpp:463
ObjCryst::ZScatterer::EndOptimization
virtual void EndOptimization()
This should be called by any optimization class at the end of an optimization.
Definition: ZScatterer.cpp:1146
ObjCryst::ZScatterer::GetComponentName
virtual string GetComponentName(const int i) const
Name for the i-th component of this scatterer.
Definition: ZScatterer.cpp:417
ObjCryst::ZScatterer::GetPsi
REAL GetPsi() const
Access to psi parameter (overall orientation of the scatterer)
Definition: ZScatterer.cpp:432
ObjCryst::ZScatterer::SetPhi
void SetPhi(const REAL)
Access to phi parameter (overall orientation of the scatterer)
Definition: ZScatterer.cpp:434
ObjCryst::ZScatterer::GetZDihedralAngleAtom
long GetZDihedralAngleAtom(const int i) const
Index of the 3rd atom used to define the i-th atom in the Z-Matrix (the one from which the dihedral a...
Definition: ZScatterer.cpp:447
ObjCryst::ZScatterer::GetZAtomY
REAL GetZAtomY(const int i) const
Get the Y fractionnal coordinate of atom i.
Definition: ZScatterer.cpp:439
ObjCryst::ZScatterer::mZAtomRegistry
ObjRegistry< ZAtom > mZAtomRegistry
Registry for ZAtoms in this Scatterer.
Definition: ZScatterer.h:441
ObjCryst::ZScatterer::GetChi
REAL GetChi() const
Access to chi parameter (overall orientation of the scatterer)
Definition: ZScatterer.cpp:431
ObjCryst::ZScatterer::GetZAngle
REAL GetZAngle(const int i) const
Const access to the angle parameter, for the i-th row in the Z-Matrix.
Definition: ZScatterer.cpp:452
ObjCryst::ZScatterer::GetGeneGroup
virtual void GetGeneGroup(const RefinableObj &obj, CrystVector_uint &groupIndex, unsigned int &firstGroup) const
Get the gene group assigned to each parameter.
Definition: ZScatterer.cpp:1103
ObjCryst::ZScatterer::mUseGlobalScattPow
bool mUseGlobalScattPow
Does the ZScatterer use a global scattering power ?
Definition: ZScatterer.h:451
ObjCryst::ZScatterer::GetZCoord
const CrystVector_REAL & GetZCoord() const
Get the list of all ZAtom cartesian x coordinates.
Definition: ZScatterer.cpp:1140
ObjCryst::ZScatterer::m3DDisplayIndex
CrystMatrix_long m3DDisplayIndex
For 3D display of the structure, bonds, triangular and quadric faces can be displayed.
Definition: ZScatterer.h:395
ObjCryst::ZScatterer::mScattCompList
ScatteringComponentList mScattCompList
The list of scattering components.
Definition: ZScatterer.h:397
ObjCryst::ZScatterer::UpdateCoordinates
void UpdateCoordinates() const
Update the atom coordinates (in real units, in Angstroems).
Definition: ZScatterer.cpp:1828
ObjCryst::ZScatterer::GetYCoord
const CrystVector_REAL & GetYCoord() const
Get the list of all ZAtom cartesian x coordinates.
Definition: ZScatterer.cpp:1135
ObjCryst::ZScatterer::mComponentIndex
CrystVector_int mComponentIndex
Index of atoms in the ScatteringComponentList.
Definition: ZScatterer.h:412
ObjCryst::ZScatterer::GetCenterAtomIndex
unsigned int GetCenterAtomIndex() const
Get the index of the central atom (around which the rotation is made)
Definition: ZScatterer.cpp:1531
ObjCryst::ZScatterer::GetClassName
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
Definition: ZScatterer.cpp:311
ObjCryst::ZScatterer::Print
void Print() const
Print a single line of information about this scatterer.
Definition: ZScatterer.cpp:422
ObjCryst::ZScatterer::begin
vector< ZAtom * >::const_iterator begin() const
low-level access to the underlying vector of ZAtom begin().
Definition: ZScatterer.cpp:1538
ObjCryst::ZScatterer::AddAtom
void AddAtom(const string &name, const ScatteringPower *pow, const long atomBond, const REAL bondLength, const long atomAngle, const REAL bondAngle, const long atomDihedral, const REAL dihedralAngle, const REAL popu=1.)
Add an atom to the Zscatterer.
Definition: ZScatterer.cpp:317
ObjCryst::ZScatterer::POVRayDescription
virtual ostream & POVRayDescription(ostream &os, const CrystalPOVRayOptions &options) const
Definition: ZScatterer.cpp:469
ObjCryst::ZScatterer::SetChi
void SetChi(const REAL)
Access to chi parameter (overall orientation of the scatterer)
Definition: ZScatterer.cpp:435
ObjCryst::ZScatterer::end
vector< ZAtom * >::const_iterator end() const
low-level access to the underlying vector of ZAtom end().
Definition: ZScatterer.cpp:1543
ObjCryst::ZScatterer::SetZBondLength
void SetZBondLength(const int i, const REAL)
Access to bondlength parameter, for the i-th row in the Z-Matrix.
Definition: ZScatterer.cpp:457
ObjCryst::ZScatterer::UpdateScattCompList
void UpdateScattCompList() const
Update the scattering component list, ie compute all atom positions from the bonds/angles/dihedral an...
Definition: ZScatterer.cpp:2014
ObjCryst::ZScatterer::mNbDummyAtom
long mNbDummyAtom
Number of "dummy" atoms in the structure.
Definition: ZScatterer.h:405
ObjCryst::ZScatterer::GetZDihedralAngle
REAL GetZDihedralAngle(const int i) const
Const access to the dihedral angle parameter, for the i-th row in the Z-Matrix.
Definition: ZScatterer.cpp:454
ObjCryst::ZScatterer::GetZBondLength
REAL GetZBondLength(const int i) const
Const access to bondlength parameter, for the i-th row in the Z-Matrix.
Definition: ZScatterer.cpp:450
ObjCryst::ZScatterer::mPhiChiPsiMatrix
CrystMatrix_REAL mPhiChiPsiMatrix
Rotation matrix for the orientation of the scatterer.
Definition: ZScatterer.h:447
ObjCryst::ZScatterer::SetPsi
void SetPsi(const REAL)
Access to psi parameter (overall orientation of the scatterer)
Definition: ZScatterer.cpp:436
ObjCryst::ZScatterer::mNbAtom
long mNbAtom
Total number of atoms in the structure.
Definition: ZScatterer.h:399
ObjCryst::ZScatterer::GetZAngleAtom
long GetZAngleAtom(const int i) const
Index of the 2nd atom used to define the i-th atom in the Z-Matrix (the one from which the angle is c...
Definition: ZScatterer.cpp:444
ObjCryst::ZScatterer::SetUseGlobalScatteringPower
virtual void SetUseGlobalScatteringPower(const bool useIt)
use a Global scattering power for this scatterer ?
Definition: ZScatterer.cpp:1072
ObjCryst::ZScatterer::GetScatteringComponentList
virtual const ScatteringComponentList & GetScatteringComponentList() const
Get the list of all scattering components for this scatterer.
Definition: ZScatterer.cpp:410
ObjCryst::ZPolyhedron
ZPolyhedron: a Scatterer to describe polyhedras such as octahedron, tetrahedron, square plane,...
Definition: ZScatterer.h:490
ObjCryst::ZPolyhedron::CreateCopy
virtual ZPolyhedron * CreateCopy() const
Definition: ZScatterer.cpp:2775
ObjCryst::ZPolyhedron::ZPolyhedron
ZPolyhedron(const RegularPolyhedraType type, Crystal &cryst, const REAL x, const REAL y, const REAL z, const string &name, const ScatteringPower *centralAtomPow, const ScatteringPower *periphAtomPow, const REAL centralPeriphDist, const REAL ligandPopu=1, const REAL phi=0., const REAL chi=0., const REAL psi=0.)
ZPolyhedron constructor.
Definition: ZScatterer.cpp:2192
ObjCryst::MonteCarloObj::Optimize
virtual void Optimize(long &nbSteps, const bool silent=false, const REAL finalcost=0, const REAL maxTime=-1)
Launch optimization (a single run) for N steps.
Definition: GlobalOptimObj.cpp:690
ObjCryst::RefParType
class of refinable parameter types.
Definition: RefinableObj.h:80
ObjCryst::RefParType::IsDescendantFromOrSameAs
bool IsDescendantFromOrSameAs(const RefParType *type) const
Returns true if the parameter is a descendant of 'type'.
Definition: RefinableObj.cpp:57
ObjCryst::RefinableObjClock::Click
void Click()
Record an event for this clock (generally, the 'time' an object has been modified,...
Definition: RefinableObj.cpp:131
ObjCryst::RefinablePar
Generic class for parameters of refinable objects.
Definition: RefinableObj.h:225
ObjCryst::RefinablePar::GetValue
REAL GetValue() const
of the parameter.
Definition: RefinableObj.cpp:338
ObjCryst::RefinablePar::CopyAttributes
void CopyAttributes(const RefinablePar &)
Copy all attributes (limits, flags, etc...) from another RefinablePar object.
Definition: RefinableObj.cpp:316
ObjCryst::RefinablePar::AssignClock
void AssignClock(RefinableObjClock &clock)
Definition: RefinableObj.cpp:664
ObjCryst::RefinablePar::GetHumanValue
const REAL & GetHumanValue() const
Current value of parameter, scaled if necessary (for angles) to a human-understandable value.
Definition: RefinableObj.cpp:398
ObjCryst::RefinablePar::Mutate
void Mutate(const REAL mutateValue)
Add the given amount to the parameter current value.
Definition: RefinableObj.cpp:445
ObjCryst::RefinablePar::GetGlobalOptimStep
REAL GetGlobalOptimStep() const
Maximum step to use during Global Optimization algorithms.
Definition: RefinableObj.cpp:596
ObjCryst::ObjRegistry
Object Registry.
Definition: RefinableObj.h:645
ObjCryst::RefinableObj
Generic Refinable Object.
Definition: RefinableObj.h:784
ObjCryst::RefinableObj::EndOptimization
virtual void EndOptimization()
This should be called by any optimization class at the end of an optimization.
Definition: RefinableObj.cpp:1721
ObjCryst::RefinableObj::AddPar
void AddPar(const RefinablePar &newRefPar)
Add a refinable parameter.
Definition: RefinableObj.cpp:1455
ObjCryst::RefinableObj::SetName
virtual void SetName(const string &name)
Name of the object.
Definition: RefinableObj.cpp:1289
ObjCryst::RefinableObj::GetPar
RefinablePar & GetPar(const long i)
Access all parameters in the order they were inputted.
Definition: RefinableObj.cpp:1391
ObjCryst::RefinableObj::GetName
virtual const string & GetName() const
Name of the object.
Definition: RefinableObj.cpp:1288
ObjCryst::RefinableObj::GetNbPar
long GetNbPar() const
Total number of refinable parameter in the object.
Definition: RefinableObj.cpp:1387
ObjCryst::RefinableObj::ResetParList
void ResetParList()
Re-init the list of refinable parameters, removing all parameters.
Definition: RefinableObj.cpp:1868
ObjCryst::RefinableObj::mOptimizationDepth
int mOptimizationDepth
Is the object being refined or optimized ? if mOptimizationDepth=0, no optimization is taking place.
Definition: RefinableObj.h:1291
ObjCryst::RefinableObj::FixAllPar
void FixAllPar()
Fix All parameters.
Definition: RefinableObj.cpp:1327
ObjCryst::RefinableObj::SetLimitsRelative
void SetLimitsRelative(const string &parName, const REAL min, const REAL max)
Change the limits for a given parameter, giving relative new limits (eg giving -.1 and +....
Definition: RefinableObj.cpp:1638
ObjCryst::RefinableObj::mName
string mName
Name for this RefinableObject. Should be unique, at least in the same scope.+.
Definition: RefinableObj.h:1263
ObjCryst::RefinableObj::GlobalOptRandomMove
virtual void GlobalOptRandomMove(const REAL mutationAmplitude, const RefParType *type=gpRefParTypeObjCryst)
Make a random move of the current configuration.
Definition: RefinableObj.cpp:1773
FormatHorizVector
Format vector as horiz array:
Definition: VFNStreamFormat.h:169
ObjCryst::WXCrystObjBasic
Abstract base class for all objects in wxCryst.
Definition: wxCryst.h:128
ObjCryst::WXZScatterer
wxCryst class for ZScatterer objects
Definition: wxZScatterer.h:48