Atom.cpp

/*  ObjCryst++ Object-Oriented Crystallographic Library
    (c) 2000-2002 Vincent Favre-Nicolin vincefn@users.sourceforge.net
        2000-2001 University of Geneva (Switzerland)
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/
/*   Atom.cpp
*  source file for the Atom scatterer
*
*/
#include <cmath>
#include <typeinfo>
#include <stdio.h> //for sprintf()
 
#include "ObjCryst/ObjCryst/Atom.h"
 
#include "ObjCryst/Quirks/VFNStreamFormat.h" //simple formatting of integers, REALs..
#include "ObjCryst/Quirks/VFNDebug.h"
 
#ifdef OBJCRYST_GL
   #ifdef __DARWIN__
      #include <OpenGL/glu.h>
   #else
      #include <GL/glu.h>
   #endif
#endif
 
#ifdef __WX__CRYST__
   #include "ObjCryst/wxCryst/wxAtom.h"
#endif
 
#include <fstream>
#include <iomanip>
 
namespace ObjCryst
{
 
//
//    ATOM : the basic atom, within the crystal
//
//includes thermic factor (betas) and x,y,z,population
//
Atom::Atom()
:mScattCompList(1),mpScattPowAtom(0)
{
   VFN_DEBUG_MESSAGE("Atom::Atom()",5)
   this->InitRefParList();
   mScattCompList(0).mpScattPow=mpScattPowAtom;
}
 
Atom::Atom( const REAL x, const REAL y, const REAL z,
            const string &name,const ScatteringPower *pow)
:mScattCompList(1),mpScattPowAtom(pow)
{
   VFN_DEBUG_MESSAGE("Atom::Atom(x,y,z,name,ScatteringPower):"<<name,5)
   this->Init(x,y,z,name,pow,1);
}
 
Atom::Atom( const REAL x, const REAL y, const REAL z, const string &name,
            const ScatteringPower *pow,const REAL popu)
:mScattCompList(1),mpScattPowAtom(pow)
{
   VFN_DEBUG_MESSAGE("Atom::Atom(x,y,z,P,B,name,ScatteringPower):"<<name,5)
   this->Init(x,y,z,name,mpScattPowAtom,popu);
}
 
Atom::Atom(const Atom &old)
:Scatterer(old),mScattCompList(1),mpScattPowAtom(old.mpScattPowAtom)
{
   VFN_DEBUG_MESSAGE("Atom::Atom(&old):/Name="<<old.mName,5)
   //:TODO: Check if this is enough (copy constructor for Scatterer ??)
   this->Init(old.mXYZ(0),old.mXYZ(1),old.mXYZ(2),
              old.mName,mpScattPowAtom,
              old.mOccupancy);
   this->GetPar(mXYZ.data()).  CopyAttributes(old.GetPar(old.mXYZ.data()));
   this->GetPar(mXYZ.data()+1).CopyAttributes(old.GetPar(old.mXYZ.data()+1));
   this->GetPar(mXYZ.data()+2).CopyAttributes(old.GetPar(old.mXYZ.data()+2));
   this->GetPar(&mOccupancy).  CopyAttributes(old.GetPar(&(old.mOccupancy)));
}
 
Atom* Atom::CreateCopy() const
{
   VFN_DEBUG_MESSAGE("Atom::CreateCopy():/Name="<<mName,10)
   return new Atom(*this);
}
 
 
Atom::~Atom()
{
   VFN_DEBUG_MESSAGE("Atom::~Atom():("<<mName<<")",5)
   if(mpScattPowAtom!=0) mpScattPowAtom->DeRegisterClient(*this);
}
 
const string& Atom::GetClassName()const
{
   const static string className="Atom";
   return className;
}
 
void Atom::operator=(const Atom &rhs)
{
   VFN_DEBUG_MESSAGE("Atom::operator=():/Name="<<rhs.mName,5)
   //:TODO: Check if this is enough (copy constructor for Scatterer ??)
   mScattCompList.Reset();
   ++mScattCompList;
 
   this->Init(rhs.mXYZ(0),rhs.mXYZ(1),rhs.mXYZ(2),
              rhs.mName,rhs.mpScattPowAtom,
              rhs.mOccupancy);
   this->GetPar(mXYZ.data()).  CopyAttributes(rhs.GetPar(rhs.mXYZ.data()));
   this->GetPar(mXYZ.data()+1).CopyAttributes(rhs.GetPar(rhs.mXYZ.data()+1));
   this->GetPar(mXYZ.data()+2).CopyAttributes(rhs.GetPar(rhs.mXYZ.data()+2));
   this->GetPar(&mOccupancy).  CopyAttributes(rhs.GetPar(&(rhs.mOccupancy)));
}
 
void Atom::Init(const REAL x, const REAL y, const REAL z,
            const string &name, const ScatteringPower *pow,
            const REAL popu)
{
   VFN_DEBUG_MESSAGE("Atom::Init():"<<name,3)
   mName=name;
   mpScattPowAtom=pow;
   mScattCompList(0).mpScattPow=mpScattPowAtom;
   mXYZ(0)=x;
   mXYZ(1)=y;
   mXYZ(2)=z;
   if(0==mpScattPowAtom)
   {//Dummy atom
      VFN_DEBUG_MESSAGE("Atom::Init()Dummy Atom:/Name="<<this->GetName(),5)
      mOccupancy=0;
      return;
   }
   mpScattPowAtom->RegisterClient(*this);
   mOccupancy=popu;
 
   if(this->GetNbPar()<4) this->InitRefParList();
 
   mClockScatterer.Click();
   VFN_DEBUG_MESSAGE("Atom::Init():End.",5)
}
 
int Atom::GetNbComponent() const {return 1;}
const ScatteringComponentList& Atom::GetScatteringComponentList() const
{
   mScattCompList(0).mX=mXYZ(0);
   mScattCompList(0).mY=mXYZ(1);
   mScattCompList(0).mZ=mXYZ(2);
   mScattCompList(0).mOccupancy=mOccupancy;
   mClockScattCompList.Click();
   return mScattCompList;
}
 
string Atom::GetComponentName(const int i) const{ return this->GetName();}
 
void Atom::Print() const
{
   VFN_DEBUG_MESSAGE("Atom::Print()",1)
   if(this->IsDummy())
   {
      cout << "Atom (DUMMY) :"
           << FormatString(this->GetName(),16)  << " at : "
           << FormatFloat(this->GetX())
           << FormatFloat(this->GetY())
           << FormatFloat(this->GetZ());
   }
   else
   {
      cout << "Atom ("
           << FormatString(mpScattPowAtom->GetSymbol(),4) << ") :"
           << FormatString(this->GetName(),16)  << " at : "
           << FormatFloat(this->GetX())
           << FormatFloat(this->GetY())
           << FormatFloat(this->GetZ())
           << ", Biso="
           << FormatFloat(mpScattPowAtom->GetBiso())
           << ", Popu=" << FormatFloat(this->GetOccupancy());
   }
   cout << endl;
}
 
REAL Atom::GetMass() const
{
   if(true==this->IsDummy()) return 0;
   return mpScattPowAtom->GetBiso();
}
 
REAL Atom::GetRadius() const
{
   if(this->IsDummy()) return 0.5;
   return mpScattPowAtom->GetRadius();
}
 
ostream& Atom::POVRayDescription(ostream &os,
                                 const CrystalPOVRayOptions &options)const
{
   if(this->IsDummy()) return os;
   if(options.mShowHydrogens==false && (mpScattPowAtom->GetForwardScatteringFactor(RAD_XRAY)<1.5)) return os;
   const REAL xMin=options.mXmin; const REAL xMax=options.mXmax;
   const REAL yMin=options.mYmin; const REAL yMax=options.mYmax;
   const REAL zMin=options.mZmin; const REAL zMax=options.mZmax;
   REAL x0,y0,z0;
   x0=mXYZ(0);
   y0=mXYZ(1);
   z0=mXYZ(2);
   const REAL aa=this->GetCrystal().GetLatticePar(0);
   const REAL bb=this->GetCrystal().GetLatticePar(1);
   const REAL cc=this->GetCrystal().GetLatticePar(2);
   CrystMatrix_REAL xyzCoords ;
   xyzCoords=this->GetCrystal().GetSpaceGroup().GetAllSymmetrics(x0,y0,z0,false,false,true);
   int nbSymmetrics=xyzCoords.rows();
   os << "// Description of Atom :" << this->GetName()<< endl;
   for(int i=0;i<nbSymmetrics;i++)
   {
      x0=xyzCoords(i,0);
      y0=xyzCoords(i,1);
      z0=xyzCoords(i,2);
      x0 = fmod((float) x0,(float)1); if(x0<0) x0+=1.;
      y0 = fmod((float) y0,(float)1); if(y0<0) y0+=1.;
      z0 = fmod((float) z0,(float)1); if(z0<0) z0+=1.;
      //Generate also translated atoms near the unit cell
      CrystMatrix_int translate(27,3);
      translate=  -1,-1,-1,
                  -1,-1, 0,
                  -1,-1, 1,
                  -1, 0,-1,
                  -1, 0, 0,
                  -1, 0, 1,
                  -1, 1,-1,
                  -1, 1, 0,
                  -1, 1, 1,
                   0,-1,-1,
                   0,-1, 0,
                   0,-1, 1,
                   0, 0,-1,
                   0, 0, 0,
                   0, 0, 1,
                   0, 1,-1,
                   0, 1, 0,
                   0, 1, 1,
                   1,-1,-1,
                   1,-1, 0,
                   1,-1, 1,
                   1, 0,-1,
                   1, 0, 0,
                   1, 0, 1,
                   1, 1,-1,
                   1, 1, 0,
                   1, 1, 1;
      for(int j=0;j<translate.rows();j++)
      {
         REAL x=x0+translate(j,0);
         REAL y=y0+translate(j,1);
         REAL z=z0+translate(j,2);
         const bool isinside=((x>=xMin) && (x<=xMax)) && ((y>=yMin) && (y<=yMax)) && ((z>=zMin) && (z<=zMax));
         REAL borderdist;
         if(isinside) borderdist=0;
         else
         {
            borderdist=0;
            if(xMin>x) borderdist+=(xMin-x)*aa*(xMin-x)*aa;
            if(yMin>y) borderdist+=(yMin-y)*bb*(yMin-y)*bb;
            if(zMin>z) borderdist+=(zMin-z)*cc*(zMin-z)*cc;
            if(xMax<x) borderdist+=(xMax-x)*aa*(xMax-x)*aa;
            if(yMax<y) borderdist+=(yMax-y)*bb*(yMax-y)*bb;
            if(zMax<z) borderdist+=(zMax-z)*cc*(zMax-z)*cc;
            borderdist=sqrt(borderdist);
         }
         REAL fout=1.0;
         if(isinside==false) fout=exp(-borderdist)*this->GetCrystal().GetDynPopCorr(this,0);
         if(fout>0.001)
         {
            this->GetCrystal().FractionalToOrthonormalCoords(x,y,z);
            os << "   ObjCrystAtom("
               <<x<<","
               <<y<<","
               <<z<<","
               <<this->GetScatteringPower().GetRadius()/3.0<<","
               <<"colour_"+this->GetScatteringPower().GetName()<<","
               <<this->GetOccupancy()<<","<<fout
               <<")"<<endl;
         }
      }
   }
   return os;
}
 
#ifdef OBJCRYST_GL
void Atom::GLInitDisplayList(const bool onlyIndependentAtoms,
                             const REAL xMin,const REAL xMax,
                             const REAL yMin,const REAL yMax,
                             const REAL zMin,const REAL zMax,
                             const bool displayEnantiomer,
                             const bool displayNames,
                             const bool hideHydrogens,
                             const REAL fadeDistance,
                             const bool fullMoleculeInLimits)const
{
   VFN_DEBUG_MESSAGE("Atom::GLInitDisplayList():"<<this->GetName(),5)
   if(this->IsDummy()) return ;
   REAL en=1;
   if(displayEnantiomer==true) en=-1;
 
      const float r=mpScattPowAtom->GetColourRGB()[0];
      const float g=mpScattPowAtom->GetColourRGB()[1];
      const float b=mpScattPowAtom->GetColourRGB()[2];
      const float f=mOccupancy;
 
      const GLfloat colour0[] = {.0, .0, .0, 0.0};
      GLfloat colourChar [] = {1.0, 1.0, 1.0, 1.0};
      if((r>0.8)&&(g>0.8)&&(b>0.8))
      {
         colourChar[0] = 0.5;
         colourChar[1] = 0.5;
         colourChar[2] = 0.5;
      }
   const REAL aa=this->GetCrystal().GetLatticePar(0);
   const REAL bb=this->GetCrystal().GetLatticePar(1);
   const REAL cc=this->GetCrystal().GetLatticePar(2);
 
   if(hideHydrogens  && (mpScattPowAtom->GetForwardScatteringFactor(RAD_XRAY)<1.5)) return;
   GLUquadricObj* pQuadric = gluNewQuadric();
   if(true==onlyIndependentAtoms)
   {
      const GLfloat colourAtom [] = {r, g, b, f};
      REAL x,y,z;
      x=mXYZ(0);
      y=mXYZ(1);
      z=mXYZ(2);
      x = fmod((REAL)x,(int)1); if(x<0) x+=1.;
      y = fmod((REAL)y,(int)1); if(y<0) y+=1.;
      z = fmod((REAL)z,(int)1); if(z<0) z+=1.;
      this->GetCrystal().FractionalToOrthonormalCoords(x,y,z);
      glPushMatrix();
         glTranslatef(x*en, y, z);
         if(displayNames)
         {
            glMaterialfv(GL_FRONT, GL_AMBIENT,   colour0);
            glMaterialfv(GL_FRONT, GL_DIFFUSE,   colour0);
            glMaterialfv(GL_FRONT, GL_SPECULAR,  colour0);
            glMaterialfv(GL_FRONT, GL_EMISSION,  colourChar);
            glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
            glRasterPos3f(0,0,0);
            crystGLPrint(this->GetName());
         }
         else
         {
            glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,   colourAtom);
            glMaterialfv(GL_FRONT, GL_SPECULAR,  colour0);
            glMaterialfv(GL_FRONT, GL_EMISSION,  colour0);
            glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
            glPolygonMode(GL_FRONT, GL_FILL);
            gluSphere(pQuadric,this->GetRadius()/3.,20,20);
         }
      glPopMatrix();
   }
   else
   {
      REAL x0,y0,z0;
      x0=mXYZ(0);
      y0=mXYZ(1);
      z0=mXYZ(2);
      CrystMatrix_REAL xyzCoords ;
      xyzCoords=this->GetCrystal().GetSpaceGroup().GetAllSymmetrics(x0,y0,z0,false,false,true);
      int nbSymmetrics=xyzCoords.rows();
      for(int i=0;i<nbSymmetrics;i++)
      {
         x0=xyzCoords(i,0);
         y0=xyzCoords(i,1);
         z0=xyzCoords(i,2);
         x0 = fmod((REAL) x0,(int)1); if(x0<0) x0+=1.;
         y0 = fmod((REAL) y0,(int)1); if(y0<0) y0+=1.;
         z0 = fmod((REAL) z0,(int)1); if(z0<0) z0+=1.;
         //Generate also translated atoms near the unit cell
         CrystMatrix_int translate(27,3);
         translate=  -1,-1,-1,
                     -1,-1, 0,
                     -1,-1, 1,
                     -1, 0,-1,
                     -1, 0, 0,
                     -1, 0, 1,
                     -1, 1,-1,
                     -1, 1, 0,
                     -1, 1, 1,
                      0,-1,-1,
                      0,-1, 0,
                      0,-1, 1,
                      0, 0,-1,
                      0, 0, 0,
                      0, 0, 1,
                      0, 1,-1,
                      0, 1, 0,
                      0, 1, 1,
                      1,-1,-1,
                      1,-1, 0,
                      1,-1, 1,
                      1, 0,-1,
                      1, 0, 0,
                      1, 0, 1,
                      1, 1,-1,
                      1, 1, 0,
                      1, 1, 1;
         for(int j=0;j<translate.rows();j++)
         {
            REAL x=x0+translate(j,0);
            REAL y=y0+translate(j,1);
            REAL z=z0+translate(j,2);
            const bool isinside=((x>=xMin) && (x<=xMax)) && ((y>=yMin) && (y<=yMax)) && ((z>=zMin) && (z<=zMax));
            REAL borderdist;
            if(isinside) borderdist=0;
            else
            {
               borderdist=0;
               if(xMin>x) borderdist+=(xMin-x)*aa*(xMin-x)*aa;
               if(yMin>y) borderdist+=(yMin-y)*bb*(yMin-y)*bb;
               if(zMin>z) borderdist+=(zMin-z)*cc*(zMin-z)*cc;
               if(xMax<x) borderdist+=(xMax-x)*aa*(xMax-x)*aa;
               if(yMax<y) borderdist+=(yMax-y)*bb*(yMax-y)*bb;
               if(zMax<z) borderdist+=(zMax-z)*cc*(zMax-z)*cc;
               borderdist=sqrt(borderdist);
            }
            REAL fout=1;
            // NB about dyn pop corr: it's not taken into account for atoms inside the view range,
            // to avoid transparency for fully occupied atoms.
            // :TODO: Maybe it should for partially occupied atoms ?
            if(isinside==false)
            {
               if ((fadeDistance == 0) || borderdist>fadeDistance) fout = 0;
               else fout*=(fadeDistance-borderdist)/fadeDistance*this->GetCrystal().GetDynPopCorr(this,0);
            }
            if(fout>0.01)
            {
               const GLfloat colourAtom [] = {r, g, b, f*fout};
               this->GetCrystal().FractionalToOrthonormalCoords(x,y,z);
               glPushMatrix();
                  glTranslatef(x*en, y, z);
                  if(displayNames)
                  {
                     if(fout>0.99)
                     {
                        glMaterialfv(GL_FRONT, GL_AMBIENT,   colour0);
                        glMaterialfv(GL_FRONT, GL_DIFFUSE,   colour0);
                        glMaterialfv(GL_FRONT, GL_SPECULAR,  colour0);
                        glMaterialfv(GL_FRONT, GL_EMISSION,  colourChar);
                        glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
                        glRasterPos3f(0,0,0);
                        crystGLPrint(this->GetName());
                     }
                  }
                  else
                  {
                     glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,   colourAtom);
                     glMaterialfv(GL_FRONT, GL_SPECULAR,  colour0);
                     glMaterialfv(GL_FRONT, GL_EMISSION,  colour0);
                     glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
                     glPolygonMode(GL_FRONT, GL_FILL);
                     gluSphere(pQuadric,this->GetRadius()/3.,20,20);
                  }
               glPopMatrix();
            }
         }
      }
   }
   gluDeleteQuadric(pQuadric);
   VFN_DEBUG_MESSAGE("Atom::GLInitDisplayList():End",5)
}
#endif  // OBJCRYST_GL
 
bool Atom::IsDummy()const { if(0==mScattCompList(0).mpScattPow) return true; return false;}
 
const ScatteringPower& Atom::GetScatteringPower()const
{ return *mpScattPowAtom;}
 
void Atom::SetScatteringPower(const ScatteringPower &p)
{
   if(mpScattPowAtom!=&p)
   {
      if(mpScattPowAtom!=0) mpScattPowAtom->DeRegisterClient(*this);
      mpScattPowAtom = &p;
      mScattCompList(0).mpScattPow=mpScattPowAtom;
      mClockScatterer.Click();
      if(mpScattPowAtom!=0) mpScattPowAtom->RegisterClient(*this);
   }
}
 
void Atom::GetGeneGroup(const RefinableObj &obj,
                                CrystVector_uint & groupIndex,
                                unsigned int &first) const
{
   //One group for all translation parameters
   unsigned int posIndex=0;
   VFN_DEBUG_MESSAGE("Atom::GetGeneGroup()",4)
   for(long i=0;i<obj.GetNbPar();i++)
      for(long j=0;j<this->GetNbPar();j++)
         if(&(obj.GetPar(i)) == &(this->GetPar(j)))
         {
            if(this->GetPar(j).GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattTransl))
            {
               if(posIndex==0) posIndex=first++;
               groupIndex(i)=posIndex;
            }
            else groupIndex(i)= first++;
         }
}
 
void Atom::InitRefParList()
{
   mClockScatterer.Click();
   VFN_DEBUG_MESSAGE("Atom::InitRefParList()",5)
   this->ResetParList();
   //:TODO: Add thermic factors
   {
      RefinablePar tmp("x",&mXYZ(0),0,1.,gpRefParTypeScattTranslX,
                        REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,1.,1.);
      tmp.AssignClock(mClockScatterer);
      this->AddPar(tmp);
   }
   {
      RefinablePar tmp("y",&mXYZ(1),0,1.,gpRefParTypeScattTranslY,
                        REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,1.,1.);
      tmp.AssignClock(mClockScatterer);
      this->AddPar(tmp);
   }
   {
      RefinablePar tmp("z",&mXYZ(2),0,1.,gpRefParTypeScattTranslZ,
                        REFPAR_DERIV_STEP_ABSOLUTE,false,false,true,true,1.,1.);
      tmp.AssignClock(mClockScatterer);
      this->AddPar(tmp);
   }
   if(false==this->IsDummy())
   {
      {
         RefinablePar tmp(this->GetName()+(string)"occup",&mOccupancy,0.01,1.,
                           gpRefParTypeScattOccup,REFPAR_DERIV_STEP_ABSOLUTE,true,true);
         tmp.AssignClock(mClockScatterer);
         tmp.SetGlobalOptimStep(.2);
         this->AddPar(tmp);
      }
   }
}
#ifdef __WX__CRYST__
WXCrystObjBasic* Atom::WXCreate(wxWindow* parent)
{
   //:TODO: Check mpWXCrystObj==0
   mpWXCrystObj=new WXAtom(parent,this);
   return mpWXCrystObj;
}
#endif
 
}//namespace