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