28 #include "ObjCryst/Quirks/VFNStreamFormat.h" 
   29 #include "ObjCryst/ObjCryst/Molecule.h" 
   30 #include "ObjCryst/RefinableObj/GlobalOptimObj.h" 
   34       #include <OpenGL/glu.h> 
   41    #include "ObjCryst/wxCryst/wxMolecule.h" 
   47 #define RIGID_BODY_STRICT_EXPERIMENTAL 
   50 #undef RESTRAINT_X2_X4_X6 
   57 XYZ::XYZ(REAL x0,REAL y0,REAL z0):x(x0),y(y0),z(z0){};
 
   88    return  sqrt( (at1.GetX()-at2.GetX())
 
   89                 *(at1.GetX()-at2.GetX())
 
   90                 +(at1.GetY()-at2.GetY())
 
   91                 *(at1.GetY()-at2.GetY())
 
   92                 +(at1.GetZ()-at2.GetZ())
 
   93                 *(at1.GetZ()-at2.GetZ()) );
 
   98    const REAL x21=at1.GetX()-at2.GetX();
 
   99    const REAL y21=at1.GetY()-at2.GetY();
 
  100    const REAL z21=at1.GetZ()-at2.GetZ();
 
  101    const REAL x23=at3.GetX()-at2.GetX();
 
  102    const REAL y23=at3.GetY()-at2.GetY();
 
  103    const REAL z23=at3.GetZ()-at2.GetZ();
 
  104    const REAL norm21_norm23= sqrt( (x21*x21+y21*y21+z21*z21)
 
  105                                   *(x23*x23+y23*y23+z23*z23)+1e-6);
 
  106    const REAL angle=(x21*x23+y21*y23+z21*z23)/norm21_norm23;
 
  107    if(angle>=1)  
return 0;
 
  108    if(angle<=-1) 
return M_PI;
 
  114    const REAL x21=at1.GetX()-at2.GetX();
 
  115    const REAL y21=at1.GetY()-at2.GetY();
 
  116    const REAL z21=at1.GetZ()-at2.GetZ();
 
  118    const REAL x34=at4.GetX()-at3.GetX();
 
  119    const REAL y34=at4.GetY()-at3.GetY();
 
  120    const REAL z34=at4.GetZ()-at3.GetZ();
 
  122    const REAL x23=at3.GetX()-at2.GetX();
 
  123    const REAL y23=at3.GetY()-at2.GetY();
 
  124    const REAL z23=at3.GetZ()-at2.GetZ();
 
  127    const REAL x123= y21*z23-z21*y23;
 
  128    const REAL y123= z21*x23-x21*z23;
 
  129    const REAL z123= x21*y23-y21*x23;
 
  132    const REAL x234= -(y23*z34-z23*y34);
 
  133    const REAL y234= -(z23*x34-x23*z34);
 
  134    const REAL z234= -(x23*y34-y23*x34);
 
  136    const REAL norm123_norm234= sqrt( (x123*x123+y123*y123+z123*z123)
 
  137                                     *(x234*x234+y234*y234+z234*z234)+1e-6);
 
  139    REAL angle=(x123*x234+y123*y234+z123*z234)/norm123_norm234;
 
  140    if(angle>= 1) angle=0;
 
  143       if(angle<=-1) angle=M_PI;
 
  144       else angle=acos(angle);
 
  146    if((x21*x234 + y21*y234 + z21*z234)<0) 
return -angle;
 
  152                               const map<
MolAtom*,set<MolAtom*> > &connect,
 
  153                               set<MolAtom*> &atomlist,
const MolAtom* finalAtom)
 
  155    const pair<set<MolAtom*>::iterator,
bool> status=atomlist.insert(atom);
 
  156    if(
false==status.second) 
return;
 
  157    if(finalAtom==atom) 
return;
 
  158    map<MolAtom*,set<MolAtom*> >::const_iterator c=connect.find(atom);
 
  159    set<MolAtom*>::const_iterator pos;
 
  160    for(pos=c->second.begin();pos!=c->second.end();++pos)
 
  167                               const map<
MolAtom*,set<MolAtom*> > &connect,
 
  168                               map<MolAtom*,unsigned long> &atomlist,
const unsigned long maxdepth,
unsigned long depth)
 
  170    if(atomlist.count(atom)>0)
 
  171      if(atomlist[atom]<=depth) 
return;
 
  172    atomlist[atom]=depth;
 
  173    if(depth==maxdepth) 
return;
 
  174    map<MolAtom*,set<MolAtom*> >::const_iterator c=connect.find(atom);
 
  175    set<MolAtom*>::const_iterator pos;
 
  176    for(pos=c->second.begin();pos!=c->second.end();++pos)
 
  188 MolAtom::MolAtom(
const REAL x, 
const REAL y, 
const REAL z,
 
  190 mName(name),mX(x),mY(y),mZ(z),mOccupancy(1.),mpScattPow(pPow),mpMol(&parent),mIsInRing(false)
 
  195    VFN_DEBUG_MESSAGE(
"MolAtom::MolAtom()",4)
 
  205 void MolAtom::SetName(
const string &name)
 
  218       catch(
const ObjCrystException &except)
 
  220          cerr<<
"MolAtom::SetName(): Atom parameters not yet declared in a Molecule ?"<<endl;
 
  225 const string& MolAtom::GetName()
const{
return mName;}
 
  226       string& MolAtom::GetName()     {
return mName;}
 
  228 const Molecule& MolAtom::GetMolecule()
const{
return *
mpMol;}
 
  229       Molecule& MolAtom::GetMolecule()     {
return *
mpMol;}
 
  231 const REAL& MolAtom::X()
const{
return mX;}
 
  232 const REAL& MolAtom::Y()
const{
return mY;}
 
  233 const REAL& MolAtom::Z()
const{
return mZ;}
 
  235 REAL& MolAtom::X(){
return mX;}
 
  236 REAL& MolAtom::Y(){
return mY;}
 
  237 REAL& MolAtom::Z(){
return mZ;}
 
  239 REAL MolAtom::GetX()
const{
return mX;}
 
  240 REAL MolAtom::GetY()
const{
return mY;}
 
  241 REAL MolAtom::GetZ()
const{
return mZ;}
 
  242 REAL MolAtom::GetOccupancy()
const{
return mOccupancy;}
 
  247 void MolAtom::SetOccupancy(
const REAL a){ 
mOccupancy=a;}
 
  251 void MolAtom::SetScatteringPower(
const ScatteringPower& pow){
mpScattPow=&pow;}
 
  253 void MolAtom::XMLOutput(ostream &os,
int indent)
const 
  255    VFN_DEBUG_ENTRY(
"MolAtom::XMLOutput()",4)
 
  256    for(
int i=0;i<indent;i++) os << "  " ;
 
  257    XMLCrystTag tag("Atom",false,true);
 
  258    tag.AddAttribute("Name",this->GetName());
 
  259    if(!this->
IsDummy())tag.AddAttribute("ScattPow",this->GetScatteringPower().GetName());
 
  262       ss.precision(os.precision());
 
  264       tag.AddAttribute(
"x",ss.str());
 
  268       ss.precision(os.precision());
 
  270       tag.AddAttribute(
"y",ss.str());
 
  274       ss.precision(os.precision());
 
  276       tag.AddAttribute(
"z",ss.str());
 
  280       ss.precision(os.precision());
 
  282       tag.AddAttribute(
"Occup",ss.str());
 
  285    VFN_DEBUG_EXIT(
"MolAtom::XMLOutput()",4)
 
  288 void MolAtom::XMLInput(istream &is,
const XMLCrystTag &tag)
 
  290    VFN_DEBUG_ENTRY(
"MolAtom::XMLInput()",7)
 
  292    for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
 
  294       if(
"Name"==tag.GetAttributeName(i))
 
  296          name=tag.GetAttributeValue(i);
 
  298       if(
"ScattPow"==tag.GetAttributeName(i))
 
  302       if(
"x"==tag.GetAttributeName(i))
 
  304          stringstream ss(tag.GetAttributeValue(i));
 
  307       if(
"y"==tag.GetAttributeName(i))
 
  309          stringstream ss(tag.GetAttributeValue(i));
 
  312       if(
"z"==tag.GetAttributeName(i))
 
  314          stringstream ss(tag.GetAttributeValue(i));
 
  317       if(
"Occup"==tag.GetAttributeName(i))
 
  319          stringstream ss(tag.GetAttributeValue(i));
 
  324    VFN_DEBUG_EXIT(
"MolAtom::XMLInput()",7)
 
  328 bool MolAtom::IsInRing()
const{
return mIsInRing;}
 
  331 WXCrystObjBasic* MolAtom::WXCreate(wxWindow* parent)
 
  333    VFN_DEBUG_ENTRY(
"MolAtom::WXCreate()",5)
 
  334    mpWXCrystObj=new WXMolAtom(parent,this);
 
  335    VFN_DEBUG_EXIT("
MolAtom::WXCreate()",5)
 
  338 WXCrystObjBasic* 
MolAtom::WXGet(){
return mpWXCrystObj;}
 
  339 void MolAtom::WXDelete(){
if(0!=mpWXCrystObj) 
delete mpWXCrystObj;mpWXCrystObj=0;}
 
  340 void MolAtom::WXNotifyDelete(){mpWXCrystObj=0;}
 
  348                  const REAL length0, 
const REAL sigma, 
const REAL delta,
 
  349                  Molecule &parent,
const REAL bondOrder):
 
  350 mAtomPair(make_pair(&atom1,&atom2)),
 
  351 mLength0(length0),mDelta(delta),mSigma(sigma),
 
  352 mBondOrder(bondOrder),mIsFreeTorsion(false),
mpMol(&parent)
 
  366       Molecule& MolBond::GetMolecule()     {
return *
mpMol;}
 
  369 {
return this->GetAtom1().GetName()+
"-"+this->GetAtom2().GetName();}
 
  371 void MolBond::XMLOutput(ostream &os,
int indent)
const 
  373    VFN_DEBUG_ENTRY(
"MolBond::XMLOutput()",4)
 
  374    for(
int i=0;i<indent;i++) os << "  " ;
 
  376    tag.AddAttribute("Atom1",mAtomPair.first->GetName());
 
  377    tag.AddAttribute("Atom2",mAtomPair.second->GetName());
 
  380       ss.precision(os.precision());
 
  382       tag.AddAttribute(
"Length",ss.str());
 
  386       ss.precision(os.precision());
 
  388       tag.AddAttribute(
"Delta",ss.str());
 
  392       ss.precision(os.precision());
 
  394       tag.AddAttribute(
"Sigma",ss.str());
 
  398       ss.precision(os.precision());
 
  400       tag.AddAttribute(
"BondOrder",ss.str());
 
  404       ss.precision(os.precision());
 
  406       tag.AddAttribute(
"FreeTorsion",ss.str());
 
  409    VFN_DEBUG_EXIT(
"MolBond::XMLOutput()",4)
 
  412 void MolBond::XMLInput(istream &is,
const XMLCrystTag &tag)
 
  414    VFN_DEBUG_ENTRY(
"MolBond::XMLInput():",7)
 
  415    for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
 
  417       if(
"Atom1"==tag.GetAttributeName(i))
 
  419          mAtomPair.first=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
 
  421       if(
"Atom2"==tag.GetAttributeName(i))
 
  423          mAtomPair.second=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
 
  425       if(
"Length"==tag.GetAttributeName(i))
 
  427          stringstream ss(tag.GetAttributeValue(i));
 
  430       if(
"Delta"==tag.GetAttributeName(i))
 
  432          stringstream ss(tag.GetAttributeValue(i));
 
  435       if(
"Sigma"==tag.GetAttributeName(i))
 
  437          stringstream ss(tag.GetAttributeValue(i));
 
  440       if(
"BondOrder"==tag.GetAttributeName(i))
 
  442          stringstream ss(tag.GetAttributeValue(i));
 
  445       if(
"FreeTorsion"==tag.GetAttributeName(i))
 
  447          stringstream ss(tag.GetAttributeValue(i));
 
  451    VFN_DEBUG_EXIT(
"MolBond::XMLInput():",7)
 
  458    if(!recalc) 
return mLLK;
 
  459    VFN_DEBUG_ENTRY(
"MolBond::GetLogLikelihood():",2)
 
  462    const REAL x=this->GetAtom2().GetX()-this->GetAtom1().GetX();
 
  463    const REAL y=this->GetAtom2().GetY()-this->GetAtom1().GetY();
 
  464    const REAL z=this->GetAtom2().GetZ()-this->GetAtom1().GetZ();
 
  465    const REAL length=sqrt(abs(x*x+y*y+z*z));
 
  469       const REAL tmp2=1/(length+1e-10);
 
  485    mLLK=length-(mLength0+mDelta);
 
  490       #ifdef RESTRAINT_X2_X4_X6 
  495       mLLK=mLLK2*(1+mLLK2);
 
  500       VFN_DEBUG_EXIT(
"MolBond::GetLogLikelihood():",2)
 
  503    mLLK=length-(mLength0-mDelta);
 
  507       #ifdef RESTRAINT_X2_X4_X6 
  508       const float mLLK2=mLLK*
mLLK;
 
  512       mLLK=mLLK2*(1+mLLK2);
 
  517       VFN_DEBUG_EXIT(
"MolBond::GetLogLikelihood():",2)
 
  530    map<const MolAtom*,XYZ>::const_iterator pos;
 
  531    pos=m.find(mAtomPair.first);
 
  536    pos=m.find(mAtomPair.second);
 
  538       d+= pos->second.x*mDerivAtom2.x
 
  539          +pos->second.y*mDerivAtom2.y
 
  540          +pos->second.z*mDerivAtom2.z;
 
  541    if(llk) 
return mDerivLLKCoeff*d;
 
  548    map<MolAtom*,XYZ>::iterator pos;
 
  549    pos=m.find(mAtomPair.first);
 
  556    pos=m.find(mAtomPair.second);
 
  565    cout<<this->
GetName()<<
" :LLK"<<endl;
 
  566    for(map<MolAtom*,XYZ>::const_iterator pos=m.begin();pos!=m.end();++pos)
 
  569       sprintf(buf,
"%10s Grad LLK= (%8.3f %8.3f %8.3f)",
 
  570             pos->first->GetName().c_str(),pos->second.x,pos->second.y,pos->second.z);
 
  576 const MolAtom& MolBond::GetAtom1()
const{
return *(mAtomPair.first);}
 
  577 const MolAtom& MolBond::GetAtom2()
const{
return *(mAtomPair.second);}
 
  578 MolAtom& MolBond::GetAtom1(){
return *(mAtomPair.first);}
 
  579 MolAtom& MolBond::GetAtom2(){
return *(mAtomPair.second);}
 
  580 void MolBond::SetAtom1(MolAtom &at){mAtomPair.first =&at;}
 
  581 void MolBond::SetAtom2(MolAtom &at){mAtomPair.second=&at;}
 
  582 REAL MolBond::GetLength()
const 
  587 REAL MolBond::GetLength0()
const{
return mLength0;}
 
  588 REAL MolBond::GetLengthDelta()
const{
return mDelta;}
 
  589 REAL MolBond::GetLengthSigma()
const{
return mSigma;}
 
  590 REAL MolBond::GetBondOrder()
const{
return mBondOrder;}
 
  592 REAL& MolBond::Length0(){
return mLength0;}
 
  593 REAL& MolBond::LengthDelta(){
return mDelta;}
 
  594 REAL& MolBond::LengthSigma(){
return mSigma;}
 
  595 REAL& MolBond::BondOrder(){
return mBondOrder;}
 
  597 void MolBond::SetLength0(
const REAL a){mLength0=a;}
 
  598 void MolBond::SetLengthDelta(
const REAL a){mDelta=a;}
 
  599 void MolBond::SetLengthSigma(
const REAL a){mSigma=a;}
 
  600 void MolBond::SetBondOrder(
const REAL a){mBondOrder=a;}
 
  602 bool MolBond::IsFreeTorsion()
const{
return mIsFreeTorsion;}
 
  603 void MolBond::SetFreeTorsion(
const bool isFreeTorsion)
 
  605    if(mIsFreeTorsion==isFreeTorsion) 
return;
 
  606    mIsFreeTorsion=isFreeTorsion;
 
  610 WXCrystObjBasic* MolBond::WXCreate(wxWindow* parent)
 
  612    VFN_DEBUG_ENTRY(
"MolBond::WXCreate()",5)
 
  613    mpWXCrystObj=new WXMolBond(parent,this);
 
  614    VFN_DEBUG_EXIT("
MolBond::WXCreate()",5)
 
  617 WXCrystObjBasic* 
MolBond::WXGet(){
return mpWXCrystObj;}
 
  618 void MolBond::WXDelete(){
if(0!=mpWXCrystObj) 
delete mpWXCrystObj;mpWXCrystObj=0;}
 
  619 void MolBond::WXNotifyDelete(){mpWXCrystObj=0;}
 
  627                            const REAL angle, 
const REAL sigma, 
const REAL delta,
 
  629 mAngle0(angle),mDelta(delta),mSigma(sigma),
mpMol(&parent)
 
  646 const Molecule& MolBondAngle::GetMolecule()
const{
return *
mpMol;}
 
  647       Molecule& MolBondAngle::GetMolecule()     {
return *
mpMol;}
 
  649 string MolBondAngle::GetName()
const 
  651    return this->GetAtom1().GetName()+
"-" 
  652          +this->GetAtom2().GetName()+
"-" 
  653          +this->GetAtom3().GetName();
 
  656 void MolBondAngle::XMLOutput(ostream &os,
int indent)
const 
  658    VFN_DEBUG_ENTRY(
"MolBondAngle::XMLOutput()",4)
 
  659    for(
int i=0;i<indent;i++) os << "  " ;
 
  660    XMLCrystTag tag("BondAngle",false,true);
 
  661    tag.AddAttribute("Atom1",this->GetAtom1().GetName());
 
  662    tag.AddAttribute("Atom2",this->GetAtom2().GetName());
 
  663    tag.AddAttribute("Atom3",this->GetAtom3().GetName());
 
  666       ss.precision(os.precision());
 
  667       ss <<mAngle0*RAD2DEG;
 
  668       tag.AddAttribute(
"Angle",ss.str());
 
  672       ss.precision(os.precision());
 
  674       tag.AddAttribute(
"Delta",ss.str());
 
  678       ss.precision(os.precision());
 
  680       tag.AddAttribute(
"Sigma",ss.str());
 
  683    VFN_DEBUG_EXIT(
"MolBondAngle::XMLOutput()",4)
 
  686 void MolBondAngle::XMLInput(istream &is,
const XMLCrystTag &tag)
 
  688    VFN_DEBUG_ENTRY(
"MolBondAngle::XMLInput():",4)
 
  690    for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
 
  692       if(
"Atom1"==tag.GetAttributeName(i))
 
  696       if(
"Atom2"==tag.GetAttributeName(i))
 
  700       if(
"Atom3"==tag.GetAttributeName(i))
 
  704       if(
"Angle"==tag.GetAttributeName(i))
 
  706          stringstream ss(tag.GetAttributeValue(i));
 
  710       if(
"Delta"==tag.GetAttributeName(i))
 
  712          stringstream ss(tag.GetAttributeValue(i));
 
  716       if(
"Sigma"==tag.GetAttributeName(i))
 
  718          stringstream ss(tag.GetAttributeValue(i));
 
  723    VFN_DEBUG_EXIT(
"MolBondAngle::XMLInput():",4)
 
  725 REAL& MolBondAngle::Angle0()
 
  729 REAL& MolBondAngle::AngleDelta(){
return mDelta;}
 
  730 REAL& MolBondAngle::AngleSigma(){
return mSigma;}
 
  732 REAL MolBondAngle::GetAngle0()
const{
return mAngle0;}
 
  733 REAL MolBondAngle::GetAngleDelta()
const{
return mDelta;}
 
  734 REAL MolBondAngle::GetAngleSigma()
const{
return mSigma;}
 
  736 void MolBondAngle::SetAngle0(
const REAL angle){mAngle0=angle;}
 
  737 void MolBondAngle::SetAngleDelta(
const REAL delta){mDelta=delta;}
 
  738 void MolBondAngle::SetAngleSigma(
const REAL sigma){mSigma=sigma;}
 
  740 REAL MolBondAngle::GetAngle()
const 
  742    return GetBondAngle(this->GetAtom1(),this->GetAtom2(),this->GetAtom3());
 
  749    if(!recalc) 
return mLLK;
 
  750    VFN_DEBUG_ENTRY(
"MolBondAngle::GetLogLikelihood():",2)
 
  753    const REAL x21=this->GetAtom1().GetX()-this->GetAtom2().GetX();
 
  754    const REAL y21=this->GetAtom1().GetY()-this->GetAtom2().GetY();
 
  755    const REAL z21=this->GetAtom1().GetZ()-this->GetAtom2().GetZ();
 
  756    const REAL x23=this->GetAtom3().GetX()-this->GetAtom2().GetX();
 
  757    const REAL y23=this->GetAtom3().GetY()-this->GetAtom2().GetY();
 
  758    const REAL z23=this->GetAtom3().GetZ()-this->GetAtom2().GetZ();
 
  760    const REAL n1=sqrt(abs(x21*x21+y21*y21+z21*z21));
 
  761    const REAL n3=sqrt(abs(x23*x23+y23*y23+z23*z23));
 
  762    const REAL p=x21*x23+y21*y23+z21*z23;
 
  764    const REAL a0=p/(n1*n3+1e-10);
 
  769       if(a0<=-1) angle=M_PI;
 
  775       const REAL s=1/(sqrt(1-a0*a0+1e-6));
 
  776       const REAL s0=-s/(n1*n3+1e-10);
 
  777       const REAL s1= s*p/(n3*n1*n1*n1+1e-10);
 
  778       const REAL s3= s*p/(n1*n3*n3*n3+1e-10);
 
  783       mDerivAtom3.x=s0*x21+s3*x23;
 
  784       mDerivAtom3.y=s0*y21+s3*y23;
 
  785       mDerivAtom3.z=s0*z21+s3*z23;
 
  799    mLLK=angle-(mAngle0+mDelta);
 
  803       #ifdef RESTRAINT_X2_X4_X6 
  804       const float mLLK2=mLLK*
mLLK;
 
  808       mLLK=mLLK2*(1+mLLK2);
 
  813       VFN_DEBUG_EXIT(
"MolBondAngle::GetLogLikelihood():",2)
 
  816    mLLK=angle-(mAngle0-mDelta);
 
  820       #ifdef RESTRAINT_X2_X4_X6 
  821       const float mLLK2=mLLK*
mLLK;
 
  825       mLLK=mLLK2*(1+mLLK2);
 
  830       VFN_DEBUG_EXIT(
"MolBondAngle::GetLogLikelihood():",2)
 
  843    map<const MolAtom*,XYZ>::const_iterator pos;
 
  851       d+= pos->second.x*mDerivAtom2.x
 
  852          +pos->second.y*mDerivAtom2.y
 
  853          +pos->second.z*mDerivAtom2.z;
 
  856       d+= pos->second.x*mDerivAtom3.x
 
  857          +pos->second.y*mDerivAtom3.y
 
  858          +pos->second.z*mDerivAtom3.z;
 
  859    if(llk) 
return mDerivLLKCoeff*d;
 
  866    map<MolAtom*,XYZ>::iterator pos;
 
  890    cout<<this->GetName()<<
" :LLK"<<endl;
 
  891    for(map<MolAtom*,XYZ>::const_iterator pos=m.begin();pos!=m.end();++pos)
 
  894       sprintf(buf,
"%10s Grad LLK= (%8.3f %8.3f %8.3f)",
 
  895             pos->first->GetName().c_str(),pos->second.x,pos->second.y,pos->second.z);
 
  901 const MolAtom& MolBondAngle::GetAtom1()
const{
return *(
mvpAtom[0]);}
 
  902 const MolAtom& MolBondAngle::GetAtom2()
const{
return *(
mvpAtom[1]);}
 
  903 const MolAtom& MolBondAngle::GetAtom3()
const{
return *(
mvpAtom[2]);}
 
  904 MolAtom& MolBondAngle::GetAtom1(){
return *(
mvpAtom[0]);}
 
  905 MolAtom& MolBondAngle::GetAtom2(){
return *(
mvpAtom[1]);}
 
  906 MolAtom& MolBondAngle::GetAtom3(){
return *(
mvpAtom[2]);}
 
  907 void MolBondAngle::SetAtom1(MolAtom& at){
mvpAtom[0]=&at;}
 
  908 void MolBondAngle::SetAtom2(MolAtom& at){
mvpAtom[1]=&at;}
 
  909 void MolBondAngle::SetAtom3(MolAtom& at){
mvpAtom[2]=&at;}
 
  914 WXCrystObjBasic* MolBondAngle::WXCreate(wxWindow* parent)
 
  916    VFN_DEBUG_ENTRY(
"MolBondAngle::WXCreate()",5)
 
  917    mpWXCrystObj=new WXMolBondAngle(parent,this);
 
  921 WXCrystObjBasic* 
MolBondAngle::WXGet(){
return mpWXCrystObj;}
 
  922 void MolBondAngle::WXDelete(){
if(0!=mpWXCrystObj) 
delete mpWXCrystObj;mpWXCrystObj=0;}
 
  923 void MolBondAngle::WXNotifyDelete(){mpWXCrystObj=0;}
 
  932                                    const REAL angle, 
const REAL sigma, 
const REAL delta,
 
  934 mAngle0(angle),mDelta(delta),mSigma(sigma),
mpMol(&parent)
 
  939    VFN_DEBUG_ENTRY(
"MolDihedralAngle::MolDihedralAngle()",5)
 
  945    mAngle0=fmod((REAL)mAngle0,(REAL)(2*M_PI));
 
  946    if(mAngle0<(-M_PI)) mAngle0+=2*M_PI;
 
  947    if(mAngle0>M_PI) mAngle0-=2*M_PI;
 
  948    VFN_DEBUG_EXIT(
"MolDihedralAngle::MolDihedralAngle()",5)
 
  958 const Molecule& MolDihedralAngle::GetMolecule()
const{
return *
mpMol;}
 
  959       Molecule& MolDihedralAngle::GetMolecule()     {
return *
mpMol;}
 
  961 string MolDihedralAngle::GetName()
const 
  963    return this->GetAtom1().GetName()+
"-" 
  964          +this->GetAtom2().GetName()+
"-" 
  965          +this->GetAtom3().GetName()+
"-" 
  966          +this->GetAtom4().GetName();
 
  969 void MolDihedralAngle::XMLOutput(ostream &os,
int indent)
const 
  971    VFN_DEBUG_ENTRY(
"MolDihedralAngle::XMLOutput()",4)
 
  972    for(
int i=0;i<indent;i++) os << "  " ;
 
  973    XMLCrystTag tag("DihedralAngle",false,true);
 
  974    tag.AddAttribute("Atom1",this->GetAtom1().GetName());
 
  975    tag.AddAttribute("Atom2",this->GetAtom2().GetName());
 
  976    tag.AddAttribute("Atom3",this->GetAtom3().GetName());
 
  977    tag.AddAttribute("Atom4",this->GetAtom4().GetName());
 
  980       ss.precision(os.precision());
 
  981       ss <<mAngle0*RAD2DEG;
 
  982       tag.AddAttribute(
"Angle",ss.str());
 
  986       ss.precision(os.precision());
 
  988       tag.AddAttribute(
"Delta",ss.str());
 
  992       ss.precision(os.precision());
 
  994       tag.AddAttribute(
"Sigma",ss.str());
 
  997    VFN_DEBUG_EXIT(
"MolDihedralAngle::XMLOutput()",4)
 
 1000 void MolDihedralAngle::XMLInput(istream &is,
const XMLCrystTag &tag)
 
 1002    VFN_DEBUG_ENTRY(
"MolDihedralAngle::XMLInput():",5)
 
 1004    for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
 
 1006       if(
"Atom1"==tag.GetAttributeName(i))
 
 1008          mvpAtom[0]=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
 
 1010       if(
"Atom2"==tag.GetAttributeName(i))
 
 1012          mvpAtom[1]=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
 
 1014       if(
"Atom3"==tag.GetAttributeName(i))
 
 1016          mvpAtom[2]=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
 
 1018       if(
"Atom4"==tag.GetAttributeName(i))
 
 1020          mvpAtom[3]=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
 
 1022       if(
"Angle"==tag.GetAttributeName(i))
 
 1024          stringstream ss(tag.GetAttributeValue(i));
 
 1028       if(
"Delta"==tag.GetAttributeName(i))
 
 1030          stringstream ss(tag.GetAttributeValue(i));
 
 1034       if(
"Sigma"==tag.GetAttributeName(i))
 
 1036          stringstream ss(tag.GetAttributeValue(i));
 
 1041    VFN_DEBUG_EXIT(
"MolDihedralAngle::XMLInput():",5)
 
 1044 REAL MolDihedralAngle::GetAngle()
const 
 1047    const REAL angle=
GetDihedralAngle(this->GetAtom1(),this->GetAtom2(),this->GetAtom3(),this->GetAtom4());
 
 1048    if((angle-mAngle0)>M_PI) 
return angle-2*M_PI;
 
 1049    else if((angle-mAngle0)<(-M_PI)) 
return angle+2*M_PI;
 
 1054 REAL& MolDihedralAngle::Angle0(){
return mAngle0;}
 
 1055 REAL& MolDihedralAngle::AngleDelta(){
return mDelta;}
 
 1056 REAL& MolDihedralAngle::AngleSigma(){
return mSigma;}
 
 1058 REAL MolDihedralAngle::GetAngle0()
const{
return mAngle0;}
 
 1059 REAL MolDihedralAngle::GetAngleDelta()
const{
return mDelta;}
 
 1060 REAL MolDihedralAngle::GetAngleSigma()
const{
return mSigma;}
 
 1062 void MolDihedralAngle::SetAngle0(
const REAL angle){mAngle0=angle;}
 
 1063 void MolDihedralAngle::SetAngleDelta(
const REAL delta){mDelta=delta;}
 
 1064 void MolDihedralAngle::SetAngleSigma(
const REAL sigma){mSigma=sigma;}
 
 1070    if(!recalc) 
return mLLK;
 
 1071    VFN_DEBUG_ENTRY(
"MolDihedralAngle::GetLogLikelihood():",2)
 
 1073    const REAL x21=this->GetAtom1().GetX()-this->GetAtom2().GetX();
 
 1074    const REAL y21=this->GetAtom1().GetY()-this->GetAtom2().GetY();
 
 1075    const REAL z21=this->GetAtom1().GetZ()-this->GetAtom2().GetZ();
 
 1077    const REAL x34=this->GetAtom4().GetX()-this->GetAtom3().GetX();
 
 1078    const REAL y34=this->GetAtom4().GetY()-this->GetAtom3().GetY();
 
 1079    const REAL z34=this->GetAtom4().GetZ()-this->GetAtom3().GetZ();
 
 1081    const REAL x23=this->GetAtom3().GetX()-this->GetAtom2().GetX();
 
 1082    const REAL y23=this->GetAtom3().GetY()-this->GetAtom2().GetY();
 
 1083    const REAL z23=this->GetAtom3().GetZ()-this->GetAtom2().GetZ();
 
 1086    const REAL x123= y21*z23-z21*y23;
 
 1087    const REAL y123= z21*x23-x21*z23;
 
 1088    const REAL z123= x21*y23-y21*x23;
 
 1091    const REAL x234= -(y23*z34-z23*y34);
 
 1092    const REAL y234= -(z23*x34-x23*z34);
 
 1093    const REAL z234= -(x23*y34-y23*x34);
 
 1095    const REAL n123= sqrt(x123*x123+y123*y123+z123*z123+1e-7);
 
 1096    const REAL n234= sqrt(x234*x234+y234*y234+z234*z234+1e-6);
 
 1098    const REAL p=x123*x234+y123*y234+z123*z234;
 
 1099    const REAL a0=p/(n123*n234+1e-10);
 
 1104       if(a0<=-1) angle=M_PI;
 
 1105       else angle=acos(a0);
 
 1108    if((x21*x34 + y21*y34 + z21*z34)<0) {angle=-angle;sgn=-1;}
 
 1113       const REAL s=sgn/(sqrt(1-a0*a0+1e-6));
 
 1114       const REAL s0=-s/(n123*n234+1e-10);
 
 1115       const REAL s1= s*p/(n234*n123*n123*n123+1e-10);
 
 1116       const REAL s3= s*p/(n123*n234*n234*n234+1e-10);
 
 1117       mDerivAtom1.x=s0*(-z23*y234+y23*z234)+s1*(-z23*y123+y23*z123);
 
 1118       mDerivAtom1.y=s0*(-x23*z234+z23*x234)+s1*(-x23*z123+z23*x123);
 
 1119       mDerivAtom1.z=s0*(-y23*x234+x23*y234)+s1*(-y23*x123+x23*y123);
 
 1121       mDerivAtom4.x=s0*(-z23*y123+y23*z123)+s3*(-z23*y234+y23*z234);
 
 1122       mDerivAtom4.y=s0*(-x23*z123+z23*x123)+s3*(-x23*z234+z23*x234);
 
 1123       mDerivAtom4.z=s0*(-y23*x123+x23*y123)+s3*(-y23*x234+x23*y234);
 
 1125       mDerivAtom2.x=s0*((z23-z21)*y234-y123*z34+(y21-y23)*z234+z123*y34)+s1*(y123*(z23-z21)+z123*(y21-y23))+s3*(-y234*z34+z234*y34);
 
 1126       mDerivAtom2.y=s0*((x23-x21)*z234-z123*x34+(z21-z23)*x234+x123*z34)+s1*(z123*(x23-x21)+x123*(z21-z23))+s3*(-z234*x34+x234*z34);
 
 1127       mDerivAtom2.z=s0*((y23-y21)*x234-x123*y34+(x21-x23)*y234+y123*x34)+s1*(x123*(y23-y21)+y123*(x21-x23))+s3*(-x234*y34+y234*x34);
 
 1129       mDerivAtom3.x=-(
mDerivAtom1.x+mDerivAtom2.x+mDerivAtom4.x);
 
 1130       mDerivAtom3.y=-(
mDerivAtom1.y+mDerivAtom2.y+mDerivAtom4.y);
 
 1131       mDerivAtom3.z=-(
mDerivAtom1.z+mDerivAtom2.z+mDerivAtom4.z);
 
 1140    mLLK=angle-(mAngle0+mDelta);
 
 1141    if(mLLK<(-M_PI)) mLLK += 2*M_PI;
 
 1142    if(mLLK>  M_PI ) mLLK -= 2*M_PI;
 
 1146       #ifdef RESTRAINT_X2_X4_X6 
 1147       const float mLLK2=mLLK*
mLLK;
 
 1151       mLLK=mLLK2*(1+mLLK2);
 
 1156       VFN_DEBUG_EXIT(
"MolDihedralAngle::GetLogLikelihood():",2)
 
 1159    mLLK=angle-(mAngle0-mDelta);
 
 1160    if(mLLK<(-M_PI)) mLLK += 2*M_PI;
 
 1161    if(mLLK>  M_PI ) mLLK -= 2*M_PI;
 
 1165       #ifdef RESTRAINT_X2_X4_X6 
 1166       const float mLLK2=mLLK*
mLLK;
 
 1170       mLLK=mLLK2*(1+mLLK2);
 
 1175       VFN_DEBUG_EXIT(
"MolDihedralAngle::GetLogLikelihood():",2)
 
 1188    map<const MolAtom*,XYZ>::const_iterator pos;
 
 1196       d+= pos->second.x*mDerivAtom2.x
 
 1197          +pos->second.y*mDerivAtom2.y
 
 1198          +pos->second.z*mDerivAtom2.z;
 
 1201       d+= pos->second.x*mDerivAtom3.x
 
 1202          +pos->second.y*mDerivAtom3.y
 
 1203          +pos->second.z*mDerivAtom3.z;
 
 1206       d+= pos->second.x*mDerivAtom4.x
 
 1207          +pos->second.y*mDerivAtom4.y
 
 1208          +pos->second.z*mDerivAtom4.z;
 
 1209    if(llk) 
return mDerivLLKCoeff*d;
 
 1216    map<MolAtom*,XYZ>::iterator pos;
 
 1247    cout<<this->GetName()<<
" :LLK"<<endl;
 
 1248    for(map<MolAtom*,XYZ>::const_iterator pos=m.begin();pos!=m.end();++pos)
 
 1251       sprintf(buf,
"%10s Grad LLK= (%8.3f %8.3f %8.3f)",
 
 1252             pos->first->GetName().c_str(),pos->second.x,pos->second.y,pos->second.z);
 
 1258 const MolAtom& MolDihedralAngle::GetAtom1()
const{
return *(
mvpAtom[0]);}
 
 1259 const MolAtom& MolDihedralAngle::GetAtom2()
const{
return *(
mvpAtom[1]);}
 
 1260 const MolAtom& MolDihedralAngle::GetAtom3()
const{
return *(
mvpAtom[2]);}
 
 1261 const MolAtom& MolDihedralAngle::GetAtom4()
const{
return *(
mvpAtom[3]);}
 
 1262 void MolDihedralAngle::SetAtom1(MolAtom& at){
mvpAtom[0]=&at;}
 
 1263 void MolDihedralAngle::SetAtom2(MolAtom& at){
mvpAtom[1]=&at;}
 
 1264 void MolDihedralAngle::SetAtom3(MolAtom& at){
mvpAtom[2]=&at;}
 
 1265 void MolDihedralAngle::SetAtom4(MolAtom& at){
mvpAtom[3]=&at;}
 
 1266 MolAtom& MolDihedralAngle::GetAtom1(){
return *(
mvpAtom[0]);}
 
 1267 MolAtom& MolDihedralAngle::GetAtom2(){
return *(
mvpAtom[1]);}
 
 1268 MolAtom& MolDihedralAngle::GetAtom3(){
return *(
mvpAtom[2]);}
 
 1269 MolAtom& MolDihedralAngle::GetAtom4(){
return *(
mvpAtom[3]);}
 
 1270 #ifdef __WX__CRYST__ 
 1271 WXCrystObjBasic* MolDihedralAngle::WXCreate(wxWindow* parent)
 
 1273    VFN_DEBUG_ENTRY(
"MolDihedralAngle::WXCreate()",5)
 
 1274    mpWXCrystObj=new WXMolDihedralAngle(parent,this);
 
 1276    return mpWXCrystObj;
 
 1279 void MolDihedralAngle::WXDelete(){
if(0!=mpWXCrystObj) 
delete mpWXCrystObj;mpWXCrystObj=0;}
 
 1280 void MolDihedralAngle::WXNotifyDelete(){mpWXCrystObj=0;}
 
 1287 string RigidGroup::GetName()
const 
 1289    set<MolAtom *>::const_iterator at=this->begin();
 
 1290    string name=(*at++)->GetName();
 
 1291    for(;at!=this->end();++at) name+=
", "+(*at)->GetName();
 
 1302 const std::list<MolAtom*>& MolRing::GetAtomList()
const 
 1305 std::list<MolAtom*>& MolRing::GetAtomList()
 
 1313 mQ0(1),mQ1(0),mQ2(0),mQ3(0),mIsUniQuaternion(true)
 
 1315    VFN_DEBUG_MESSAGE(
"Quaternion::Quaternion()",5)
 
 1323 mQ0(q0),mQ1(q1),mQ2(q2),mQ3(q3),mIsUniQuaternion(unit)
 
 1325    VFN_DEBUG_MESSAGE(
"Quaternion::Quaternion()",5)
 
 1329 Quaternion::~Quaternion()
 
 1331    VFN_DEBUG_MESSAGE(
"Quaternion::~Quaternion()",5)
 
 1339    VFN_DEBUG_MESSAGE(
"Quaternion::RotationQuaternion()",4)
 
 1340    const REAL s=sin(ang/2.)/sqrt(v1*v1+v2*v2+v3*v3+1e-7);
 
 1341    return Quaternion(cos(ang/2.),s*v1,s*v2,s*v3,
 
 1353       (this->Q0()*q.Q0()-this->Q1()*q.Q1()-this->Q2()*q.Q2()-this->Q3()*q.Q3(),
 
 1354        this->Q0()*q.Q1()+this->Q1()*q.Q0()+this->Q2()*q.Q3()-this->Q3()*q.Q2(),
 
 1355        this->Q0()*q.Q2()-this->Q1()*q.Q3()+this->Q2()*q.Q0()+this->Q3()*q.Q1(),
 
 1356        this->Q0()*q.Q3()+this->Q1()*q.Q2()-this->Q2()*q.Q1()+this->Q3()*q.Q0(),
false);
 
 1359 void Quaternion::operator*=(
const Quaternion &q)
 
 1363    const REAL q1=this->Q0()*q.Q1()+this->Q1()*q.Q0()+this->Q2()*q.Q3()-this->Q3()*q.Q2();
 
 1364    const REAL q2=this->Q0()*q.Q2()+this->Q2()*q.Q0()-this->Q1()*q.Q3()+this->Q3()*q.Q1();
 
 1365    const REAL q3=this->Q0()*q.Q3()+this->Q3()*q.Q0()+this->Q1()*q.Q2()-this->Q2()*q.Q1();
 
 1366    this->Q0()=   this->Q0()*q.Q0()-this->Q1()*q.Q1()-this->Q2()*q.Q2()-this->Q3()*q.Q3();
 
 1374 void Quaternion::XMLOutput(ostream &os,
int indent)
const 
 1376    VFN_DEBUG_ENTRY(
"Quaternion::XMLOutput()",4)
 
 1377    for(
int i=0;i<indent;i++) os << "  " ;
 
 1382       ss.precision(os.precision());
 
 1384       tag.AddAttribute(
"Q0",ss.str());
 
 1388       ss.precision(os.precision());
 
 1390       tag.AddAttribute(
"Q1",ss.str());
 
 1394       ss.precision(os.precision());
 
 1396       tag.AddAttribute(
"Q2",ss.str());
 
 1400       ss.precision(os.precision());
 
 1402       tag.AddAttribute(
"Q3",ss.str());
 
 1406       ss.precision(os.precision());
 
 1407       ss <<mIsUniQuaternion;
 
 1408       tag.AddAttribute(
"IsUnitQuaternion",ss.str());
 
 1411    VFN_DEBUG_EXIT(
"Quaternion::XMLOutput()",4)
 
 1414 void Quaternion::XMLInput(istream &is,
const XMLCrystTag &tag)
 
 1416    VFN_DEBUG_ENTRY(
"Quaternion::XMLInput()",5)
 
 1417    for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
 
 1419       if(
"Q0"==tag.GetAttributeName(i))
 
 1421          stringstream ss(tag.GetAttributeValue(i));
 
 1424       if(
"Q1"==tag.GetAttributeName(i))
 
 1426          stringstream ss(tag.GetAttributeValue(i));
 
 1429       if(
"Q2"==tag.GetAttributeName(i))
 
 1431          stringstream ss(tag.GetAttributeValue(i));
 
 1434       if(
"Q3"==tag.GetAttributeName(i))
 
 1436          stringstream ss(tag.GetAttributeValue(i));
 
 1439       if(
"IsUnitQuaternion"==tag.GetAttributeName(i))
 
 1441          stringstream ss(tag.GetAttributeValue(i));
 
 1442          ss >>mIsUniQuaternion;
 
 1446    VFN_DEBUG_EXIT(
"Quaternion::XMLInput()",5)
 
 1453    Quaternion P(0,v1,v2,v3,
false);
 
 1462    const REAL p0=-mQ1*v1          - mQ2*v2 - mQ3*v3;
 
 1463    const REAL p1= mQ0*v1          + mQ2*v3 - mQ3*v2;
 
 1464    const REAL p2= mQ0*v2          - mQ1*v3 + mQ3*v1;
 
 1465    const REAL p3= mQ0*v3          + mQ1*v2 - mQ2*v1;
 
 1467    v1           =-p0*mQ1 + p1*mQ0 - p2*mQ3 + p3*mQ2;
 
 1468    v2           =-p0*mQ2 + p2*mQ0 + p1*mQ3 - p3*mQ1;
 
 1469    v3           =-p0*mQ3 + p3*mQ0 - p1*mQ2 + p2*mQ1;
 
 1474    const REAL norm=sqrt( this->Q0()*this->Q0()
 
 1475                         +this->Q1()*this->Q1()
 
 1476                         +this->Q2()*this->Q2()
 
 1477                         +this->Q3()*this->Q3());
 
 1483 REAL Quaternion::GetNorm()
const 
 1484 {
return sqrt( this->Q0()*this->Q0()
 
 1485              +this->Q1()*this->Q1()
 
 1486              +this->Q2()*this->Q2()
 
 1487              +this->Q3()*this->Q3());}
 
 1489 const REAL& Quaternion::Q0()
const{
return mQ0;}
 
 1490 const REAL& Quaternion::Q1()
const{
return mQ1;}
 
 1491 const REAL& Quaternion::Q2()
const{
return mQ2;}
 
 1492 const REAL& Quaternion::Q3()
const{
return mQ3;}
 
 1493 REAL& Quaternion::Q0(){
return mQ0;}
 
 1494 REAL& Quaternion::Q1(){
return mQ1;}
 
 1495 REAL& Quaternion::Q2(){
return mQ2;}
 
 1496 REAL& Quaternion::Q3(){
return mQ3;}
 
 1502 StretchMode::~StretchMode(){}
 
 1506 mpAtom0(&at0),mpAtom1(&at1),mpBond(pBond)
 
 1509    mpMol = &(at1.GetMolecule());
 
 1512 StretchModeBondLength::~StretchModeBondLength(){}
 
 1522       const REAL n=sqrt(dx*dx+dy*dy+dz*dz+1e-7);
 
 1539       for(map<const MolBond*,REAL>::const_iterator pos=this->
mvpBrokenBond.begin();
 
 1541       for(map<const MolBondAngle*,REAL>::const_iterator pos=this->
mvpBrokenBondAngle.begin();
 
 1553       cout<<
", Translated Atoms:";
 
 1557          cout<<(*atom)->GetName()<<
",";
 
 1563                                     const bool keepCenter)
 
 1568    const REAL l=sqrt(dx*dx+dy*dy+dz*dz+1e-7);
 
 1570    const REAL change=amplitude/l;
 
 1578                                           const bool keepCenter)
 
 1585 mpAtom0(&at0),mpAtom1(&at1),mpAtom2(&at2),mpBondAngle(pBondAngle)
 
 1588    mpMol = &(at1.GetMolecule());
 
 1591 StretchModeBondAngle::~StretchModeBondAngle(){}
 
 1597       const REAL x1=
mpAtom1->GetX(),
 
 1601       const REAL dx10=
mpAtom0->GetX()-x1,
 
 1608       REAL vx=dy10*dz12-dz10*dy12,
 
 1609            vy=dz10*dx12-dx10*dz12,
 
 1610            vz=dx10*dy12-dy10*dx12;
 
 1612       const REAL n=sqrt(vx*vx+vy*vy+vz*vz+1e-10);
 
 1627          const REAL x=(*pos)->GetX()-x1,
 
 1628                     y=(*pos)->GetY()-y1,
 
 1629                     z=(*pos)->GetZ()-z1;
 
 1638       for(map<const MolBond*,REAL>::const_iterator pos=this->
mvpBrokenBond.begin();
 
 1640       for(map<const MolBondAngle*,REAL>::const_iterator pos=this->
mvpBrokenBondAngle.begin();
 
 1652       cout<<
", Rotated Atoms:";
 
 1656          cout<<(*atom)->GetName()<<
",";
 
 1662                                    const bool keepCenter)
 
 1671    const REAL vx=dy10*dz12-dz10*dy12;
 
 1672    const REAL vy=dz10*dx12-dx10*dz12;
 
 1673    const REAL vz=dx10*dy12-dy10*dx12;
 
 1678                                          const bool keepCenter)
 
 1685 mpAtom1(&at1),mpAtom2(&at2),mpDihedralAngle(pAngle)
 
 1688    mpMol = &(at1.GetMolecule());
 
 1691 StretchModeTorsion::~StretchModeTorsion(){}
 
 1698       const REAL x1=
mpAtom1->GetX(),
 
 1706       const REAL n=sqrt(vx*vx+vy*vy+vz*vz+1e-10);
 
 1715          const REAL x=(*pos)->GetX()-x1,
 
 1716                     y=(*pos)->GetY()-y1,
 
 1717                     z=(*pos)->GetZ()-z1;
 
 1726       for(map<const MolBond*,REAL>::const_iterator pos=this->
mvpBrokenBond.begin();
 
 1728       for(map<const MolBondAngle*,REAL>::const_iterator pos=this->
mvpBrokenBondAngle.begin();
 
 1740       cout<<
", Rotated Atoms:";
 
 1744          cout<<(*atom)->GetName()<<
",";
 
 1755                                        const bool keepCenter)
 
 1767 mpAtom1(&at1),mpAtom2(&at2)
 
 1770    mpMol = &(at1.GetMolecule());
 
 1773 StretchModeTwist::~StretchModeTwist(){}
 
 1779       const REAL x1=
mpAtom1->GetX(),
 
 1787       const REAL n=sqrt(vx*vx+vy*vy+vz*vz+1e-10);
 
 1796          const REAL x=(*pos)->GetX()-x1,
 
 1797                     y=(*pos)->GetY()-y1,
 
 1798                     z=(*pos)->GetZ()-z1;
 
 1807       for(map<const MolBond*,REAL>::const_iterator pos=this->
mvpBrokenBond.begin();
 
 1809       for(map<const MolBondAngle*,REAL>::const_iterator pos=this->
mvpBrokenBondAngle.begin();
 
 1821       os<<
", Rotated Atoms:";
 
 1825          os<<(*atom)->GetName()<<
",";
 
 1836                                      const bool keepCenter)
 
 1841    if((abs(dx)+abs(dy)+abs(dz))<1e-6) 
return;
 
 1842    const REAL change=(REAL)(2.*rand()-RAND_MAX)/(REAL)RAND_MAX*
mBaseAmplitude*amplitude;
 
 1855                          set<MolBondAngle*> &va,
 
 1856                          set<MolDihedralAngle*> &vd):
 
 1861    for(set<MolBond*>::iterator pos=vb.begin();pos!=vb.end();++pos)
 
 1862       mvpBond.push_back(*pos);
 
 1863    for(set<MolBondAngle*>::iterator pos=va.begin();pos!=va.end();++pos)
 
 1864       mvpBondAngle.push_back(*pos);
 
 1865    for(set<MolDihedralAngle*>::iterator pos=vd.begin();pos!=vd.end();++pos)
 
 1866       mvpDihedralAngle.push_back(*pos);
 
 1871    if(full) os<<
"MDAtomGroup: ";
 
 1872    for(set<MolAtom*>::const_iterator pos=mvpAtom.begin();pos!=mvpAtom.end();++pos)
 
 1873       os<<(*pos)->GetName()<<
" ";
 
 1876       os<<endl<<
"  Involving bond restraints:";
 
 1877       for(vector<MolBond*>::const_iterator pos=mvpBond.begin();pos!=mvpBond.end();++pos)
 
 1878          os<<(*pos)->GetName()<<
"  ";
 
 1879       os<<endl<<
"  Involving bond angle restraints:";
 
 1880       for(vector<MolBondAngle*>::const_iterator pos=mvpBondAngle.begin();pos!=mvpBondAngle.end();++pos)
 
 1881          os<<(*pos)->GetName()<<
"  ";
 
 1882       os<<endl<<
"  Involving dihedral angle restraints:";
 
 1883       for(vector<MolDihedralAngle*>::const_iterator pos=mvpDihedralAngle.begin();pos!=mvpDihedralAngle.end();++pos)
 
 1884          os<<(*pos)->GetName()<<
"  ";
 
 1895 mDeleteSubObjInDestructor(1), mBaseRotationAmplitude(M_PI*0.02), mIsSelfOptimizing(false),
 
 1896 mpCenterAtom(0), mMDMoveFreq(0.0), mMDMoveEnergy(40.), mLogLikelihoodScale(1.0)
 
 1898    VFN_DEBUG_MESSAGE(
"Molecule::Molecule()",5)
 
 1902                         gpRefParTypeScattTranslX,
 
 1903                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
 
 1910                         gpRefParTypeScattTranslY,
 
 1911                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
 
 1918                         gpRefParTypeScattTranslZ,
 
 1919                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
 
 1926                         gpRefParTypeScattOccup,
 
 1927                         REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.,1.);
 
 1934                         gpRefParTypeScattOrient,
 
 1935                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 1943                         gpRefParTypeScattOrient,
 
 1944                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 1952                         gpRefParTypeScattOrient,
 
 1953                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 1961                         gpRefParTypeScattOrient,
 
 1962                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 1969    mLocalParamSet=this->
CreateParamSet(
"saved parameters for local minimization");
 
 1983 mDeleteSubObjInDestructor(old.mDeleteSubObjInDestructor), mIsSelfOptimizing(false), mpCenterAtom(0)
 
 1985    VFN_DEBUG_ENTRY(
"Molecule::Molecule(old&)",5)
 
 1990                         gpRefParTypeScattTranslX,
 
 1991                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
 
 1998                         gpRefParTypeScattTranslY,
 
 1999                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
 
 2006                         gpRefParTypeScattTranslZ,
 
 2007                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
 
 2014                         gpRefParTypeScattOccup,
 
 2015                         REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.,1.);
 
 2022                         gpRefParTypeScattOrient,
 
 2023                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 2031                         gpRefParTypeScattOrient,
 
 2032                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 2040                         gpRefParTypeScattOrient,
 
 2041                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 2049                         gpRefParTypeScattOrient,
 
 2050                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 2056    mLocalParamSet=this->
CreateParamSet(
"saved parameters for local minimization");
 
 2072    VFN_DEBUG_EXIT(
"Molecule::Molecule(old&)",5)
 
 2077    VFN_DEBUG_ENTRY(
"Molecule::~Molecule()",5)
 
 2081           vector<MolAtom*>::iterator pos;
 
 2085           vector<MolBond*>::iterator pos;
 
 2089           vector<MolBondAngle*>::iterator pos;
 
 2093           vector<MolDihedralAngle*>::iterator pos;
 
 2097    VFN_DEBUG_EXIT(
"Molecule::~Molecule()",5)
 
 2102    VFN_DEBUG_MESSAGE(
"Molecule::CreateCopy()",5)
 
 2108    static const string className=
"Molecule";
 
 2114    if(
mName==name) 
return;
 
 2130       cerr<<
"Molecule::SetName(): parameters not yet declared in a Molecule ?"<<endl;
 
 2137    std::map<std::string,float> velts;
 
 2138    for(std::vector<MolAtom*>::const_iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 2140       if((*pos)->IsDummy()) 
continue;
 
 2142       if((*pos)->GetScatteringPower().GetClassName()==
"ScatteringPowerAtom")
 
 2143          p=dynamic_cast<const ScatteringPowerAtom*>(&((*pos)->GetScatteringPower()))->GetSymbol();
 
 2144       else p=(*pos)->GetScatteringPower().GetName();
 
 2145       if(velts.count(p)==0)
 
 2146          velts[(*pos)->GetScatteringPower().GetName()]=(*pos)->GetOccupancy();
 
 2147       else velts[(*pos)->GetScatteringPower().GetName()]+=(*pos)->GetOccupancy();
 
 2150    s<<std::setprecision(2);
 
 2151    for(std::map<std::string,float>::const_iterator pos=velts.begin();pos!=velts.end();++pos)
 
 2153       if(pos!=velts.begin()) s<<
" ";
 
 2154       float nb=pos->second;
 
 2155       if((abs(nb)-nb)<0.01) s<<pos->first<<int(round(nb));
 
 2156       else s<<pos->first<<nb;
 
 2163    VFN_DEBUG_MESSAGE(
"Molecule::Print()",5)
 
 2169    VFN_DEBUG_ENTRY(
"Molecule::XMLOutput()",4)
 
 2174    for(
int i=0;i<indent;i++) os << 
"  " ;
 
 2176    tag.AddAttribute(
"Name",
mName);
 
 2179    tag.AddAttribute(
"MDMoveFreq",sst.str());
 
 2182    tag.AddAttribute(
"MDMoveEnergy",sst.str());
 
 2185    tag.AddAttribute(
"LogLikelihoodScale",sst.str());
 
 2190    mQuat.XMLOutput(os,indent);
 
 2212       vector<MolAtom*>::const_iterator pos;
 
 2214          (*pos)->XMLOutput(os,indent);
 
 2217       vector<MolBond*>::const_iterator pos;
 
 2219          (*pos)->XMLOutput(os,indent);
 
 2222       vector<MolBondAngle*>::const_iterator pos;
 
 2224          (*pos)->XMLOutput(os,indent);
 
 2227       vector<MolDihedralAngle*>::const_iterator pos;
 
 2229          (*pos)->XMLOutput(os,indent);
 
 2232       vector<RigidGroup*>::const_iterator pos;
 
 2236          for(set<MolAtom *>::const_iterator at=(*pos)->begin();at!=(*pos)->end();++at)
 
 2237             tagg.AddAttribute(
"Atom",(*at)->GetName());
 
 2247          for(
int i=0;i<indent;i++) os << 
"  " ;
 
 2256       for(
int i=0;i<indent;i++) os << 
"  " ;
 
 2261    tag.SetIsEndTag(
true);
 
 2262    for(
int i=0;i<indent;i++) os << 
"  " ;
 
 2264    VFN_DEBUG_EXIT(
"Molecule::XMLOutput()",4)
 
 2269    VFN_DEBUG_ENTRY(
"Molecule::XMLInput()",5)
 
 2270    for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
 
 2272       if(
"Name"==tag.GetAttributeName(i))
 
 2274          mName=tag.GetAttributeValue(i);
 
 2276       if(
"MDMoveFreq"==tag.GetAttributeName(i))
 
 2280       if(
"MDMoveEnergy"==tag.GetAttributeName(i))
 
 2284       if(
"LogLikelihoodScale"==tag.GetAttributeName(i))
 
 2292       if((
"Molecule"==tagg.GetName())&&tagg.IsEndTag())
 
 2296          VFN_DEBUG_EXIT(
"Molecule::XMLInput():"<<this->
GetName(),5)
 
 2299       if(
"Quaternion"==tagg.GetName())
 
 2301          mQuat.XMLInput(is,tagg);
 
 2303       if(
"Atom"==tagg.GetName())
 
 2306          mvpAtom.back()->XMLInput(is,tagg);
 
 2308       if(
"Bond"==tagg.GetName())
 
 2310          this->
AddBond(this->GetAtom(0),this->GetAtom(1),1.5,.01,.05,1.,
false);
 
 2311          mvpBond.back()->XMLInput(is,tagg);
 
 2313       if(
"BondAngle"==tagg.GetName())
 
 2315          this->
AddBondAngle(this->GetAtom(0),this->GetAtom(1),this->GetAtom(2),1.5,.01,.05,
false);
 
 2318       if(
"DihedralAngle"==tagg.GetName())
 
 2321                                 this->GetAtom(2),this->GetAtom(3),1.5,.01,.05,
false);
 
 2324       if(
"RigidGroup"==tagg.GetName())
 
 2327          for(
unsigned int i=0;i<tagg.GetNbAttribute();i++)
 
 2328             if(
"Atom"==tagg.GetAttributeName(i))
 
 2329                s.insert(&(this->GetAtom(tagg.GetAttributeValue(i))));
 
 2332       if(
"CenterAtom"==tagg.GetName())
 
 2335          for(
unsigned int i=0;i<tagg.GetNbAttribute();i++)
 
 2336             if(
"Name"==tagg.GetAttributeName(i))
 
 2337                this->
SetCenterAtom(this->GetAtom(tagg.GetAttributeValue(i)));
 
 2339       if(
"Option"==tagg.GetName())
 
 2341          for(
unsigned int i=0;i<tagg.GetNbAttribute();i++)
 
 2342             if(
"Name"==tagg.GetAttributeName(i))
 
 2346       if(
"Par"==tagg.GetName())
 
 2348          for(
unsigned int i=0;i<tagg.GetNbAttribute();i++)
 
 2350             if(
"Name"==tagg.GetAttributeName(i))
 
 2352                if(
"x"==tagg.GetAttributeValue(i))
 
 2357                if(
"y"==tagg.GetAttributeValue(i))
 
 2362                if(
"z"==tagg.GetAttributeValue(i))
 
 2367                if(
"Occup"==tagg.GetAttributeValue(i))
 
 2376    VFN_DEBUG_EXIT(
"Molecule::XMLInput()",5)
 
 2393    TAU_PROFILE(
"Molecule::BeginOptimization()",
"void (bool,bool)",TAU_DEFAULT);
 
 2394    #if 1 // Is doing this automatically too dangerous ? 
 2395    if(  (!mIsSelfOptimizing)
 
 2400          (*fpObjCrystInformUser)(
"Optimizing initial conformation of Molecule:"+this->
GetName());
 
 2402          (*fpObjCrystInformUser)(
"");
 
 2409    clockConf=mClockAtomList;
 
 2410    if(clockConf<mClockBondList) clockConf=mClockBondList;
 
 2411    if(clockConf<mClockBondAngleList) clockConf=mClockBondAngleList;
 
 2412    if(clockConf<mClockDihedralAngleList) clockConf=mClockDihedralAngleList;
 
 2413    if(clockConf<mClockRigidGroup) clockConf=mClockRigidGroup;
 
 2414    if(clockConf<mClockAtomScattPow) clockConf=mClockAtomScattPow;
 
 2416    clockMode=mClockConnectivityTable;
 
 2417    if(clockMode<mClockRingList) clockMode=mClockRingList;
 
 2418    if(clockMode<mClockRotorGroup) clockMode=mClockRotorGroup;
 
 2419    if(clockMode<mClockFlipGroup) clockMode=mClockFlipGroup;
 
 2420    if(clockMode<mClockStretchModeBondLength) clockMode=mClockStretchModeBondLength;
 
 2421    if(clockMode<mClockStretchModeBondAngle) clockMode=mClockStretchModeBondAngle;
 
 2422    if(clockMode<mClockStretchModeTorsion) clockMode=mClockStretchModeTorsion;
 
 2423    if(clockMode<mClockStretchModeTwist) clockMode=mClockStretchModeTwist;
 
 2424    if(clockMode<mClockMDAtomGroup) clockMode=mClockMDAtomGroup;
 
 2427    if( (!mIsSelfOptimizing) && (clockMode<clockConf))
 
 2461       for(vector<MolAtom*>::iterator at=this->GetAtomList().begin();at!=this->GetAtomList().end();++at)
 
 2463          this->
GetPar(&((*at)->X())).SetIsFixed(
true);
 
 2464          this->
GetPar(&((*at)->Y())).SetIsFixed(
true);
 
 2465          this->
GetPar(&((*at)->Z())).SetIsFixed(
true);
 
 2470       for(vector<MolAtom*>::iterator at=this->GetAtomList().begin();at!=this->GetAtomList().end();++at)
 
 2472          this->
GetPar(&((*at)->X())).SetIsFixed(
false);
 
 2473          this->
GetPar(&((*at)->Y())).SetIsFixed(
false);
 
 2474          this->
GetPar(&((*at)->Z())).SetIsFixed(
false);
 
 2477    #ifdef RIGID_BODY_STRICT_EXPERIMENTAL 
 2486       (*pos)->mQuat.Q0()=1;
 
 2487       (*pos)->mQuat.Q1()=0;
 
 2488       (*pos)->mQuat.Q2()=0;
 
 2489       (*pos)->mQuat.Q3()=0;
 
 2490       (*pos)->mvIdx.clear();
 
 2491       for(set<MolAtom *>::iterator at=(*pos)->begin();at!=(*pos)->end();++at)
 
 2493          this->
GetPar(&((*at)->X())).SetIsFixed(
true);
 
 2494          this->
GetPar(&((*at)->Y())).SetIsFixed(
true);
 
 2495          this->
GetPar(&((*at)->Z())).SetIsFixed(
true);
 
 2497             if(&(this->GetAtom(i))==*at)
 
 2499               (*pos)->mvIdx.insert(i);
 
 2507    mRandomConformChangeNbTest=0;
 
 2508    mRandomConformChangeNbAccept=0;
 
 2509    mRandomConformChangeTemp=1.;
 
 2521    #ifdef RIGID_BODY_STRICT_EXPERIMENTAL 
 2525       for(set<MolAtom *>::iterator at=(*pos)->begin();at!=(*pos)->end();++at)
 
 2527          this->
GetPar(&((*at)->X())).SetIsFixed(
false);
 
 2528          this->
GetPar(&((*at)->Y())).SetIsFixed(
false);
 
 2529          this->
GetPar(&((*at)->Z())).SetIsFixed(
false);
 
 2541    TAU_PROFILE(
"Molecule::RandomizeConfiguration()",
"void ()",TAU_DEFAULT);
 
 2542    VFN_DEBUG_ENTRY(
"Molecule::RandomizeConfiguration()",4)
 
 2544    if(  (!mIsSelfOptimizing)
 
 2548       (*fpObjCrystInformUser)(
"Optimizing initial conformation of Molecule:"+this->
GetName());
 
 2550       (*fpObjCrystInformUser)(
"");
 
 2576          const REAL angle=(REAL)rand()*2.*M_PI/(REAL)RAND_MAX;
 
 2578                                pos->mvRotatedAtomList,angle);
 
 2584          const REAL angle=(REAL)rand()*2.*M_PI/(REAL)RAND_MAX;
 
 2586                                pos->mvRotatedAtomList,angle);
 
 2592          const REAL angle=(REAL)rand()*2.*M_PI/(REAL)RAND_MAX;
 
 2594                                pos->mvRotatedAtomList,angle);
 
 2597    for(list<StretchModeTorsion>::const_iterator
 
 2601       const REAL amp=2*M_PI*rand()/(REAL)RAND_MAX;
 
 2608       map<MolAtom*,XYZ> v0;
 
 2609       for(vector<MolAtom*>::iterator at=this->GetAtomList().begin();at!=this->GetAtomList().end();++at)
 
 2610          v0[*at]=
XYZ(rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5);
 
 2613                                      +this->GetBondAngleList().size()
 
 2614                                      +this->GetDihedralAngleList().size());
 
 2615       map<RigidGroup*,std::pair<XYZ,XYZ> > vr;
 
 2617                                     this->GetBondList(),
 
 2618                                     this->GetBondAngleList(),
 
 2619                                     this->GetDihedralAngleList(),
 
 2626    #ifdef RIGID_BODY_STRICT_EXPERIMENTAL 
 2634       (*pos)->mQuat.Q0()=1;
 
 2635       (*pos)->mQuat.Q1()=0;
 
 2636       (*pos)->mQuat.Q2()=0;
 
 2637       (*pos)->mQuat.Q3()=0;
 
 2642       const REAL amp=M_PI/RAND_MAX;
 
 2644                   ((2.*(REAL)rand()-(REAL)RAND_MAX)*amp,
 
 2645                    (REAL)rand(),(REAL)rand(),(REAL)rand());
 
 2647       mClockOrientation.
Click();
 
 2649    VFN_DEBUG_EXIT(
"Molecule::RandomizeConfiguration()",4)
 
 2656    if(mIsSelfOptimizing)
 
 2660       VFN_DEBUG_EXIT(
"Molecule::GlobalOptRandomMove()",4)
 
 2663    TAU_PROFILE(
"Molecule::GlobalOptRandomMove()",
"void (REAL,RefParType*)",TAU_DEFAULT);
 
 2664    TAU_PROFILE_TIMER(timer1,
"Molecule::GlobalOptRandomMove 1",
"", TAU_FIELD);
 
 2665    TAU_PROFILE_TIMER(timer2,
"Molecule::GlobalOptRandomMove 2",
"", TAU_FIELD);
 
 2666    TAU_PROFILE_TIMER(timer3,
"Molecule::GlobalOptRandomMove 3",
"", TAU_FIELD);
 
 2667    TAU_PROFILE_TIMER(timer4,
"Molecule::GlobalOptRandomMove 4",
"", TAU_FIELD);
 
 2668    TAU_PROFILE_TIMER(timer5,
"Molecule::GlobalOptRandomMove 5",
"", TAU_FIELD);
 
 2669    VFN_DEBUG_ENTRY(
"Molecule::GlobalOptRandomMove()",4)
 
 2676       &&(mFlipModel.GetChoice()==0)
 
 2679       &&(((rand()%100)==0)))
 
 2684       const unsigned long i=rand() % 
mvFlipGroup.size();
 
 2685       list<FlipGroup>::iterator pos=
mvFlipGroup.begin();
 
 2686       for(
unsigned long j=0;j<i;++j)++pos;
 
 2689       static unsigned long ctflip1=0,ctflip2=0;
 
 2690       if(pos->mvRotatedChainList.begin()->first==pos->mpAtom0)
 
 2692          cout <<
"TRYING: Flip group from atom " 
 2693               <<pos->mpAtom0->GetName()<<
",exchanging bonds with " 
 2694               <<pos->mpAtom1->GetName()<<
" and " 
 2695               <<pos->mpAtom2->GetName()<<
", resulting in a 180deg rotation of atoms : ";
 
 2696          for(set<MolAtom*>::iterator pos1=pos->mvRotatedChainList.begin()->second.begin();
 
 2697              pos1!=pos->mvRotatedChainList.begin()->second.end();++pos1)
 
 2698             cout<<(*pos1)->GetName()<<
"  ";
 
 2702          cout <<
"TRYING: Flip group with respect to: " 
 2703               <<pos->mpAtom1->GetName()<<
"-" 
 2704               <<pos->mpAtom0->GetName()<<
"-" 
 2705               <<pos->mpAtom2->GetName()<<
" : ";
 
 2706          for(list<pair<
const MolAtom *,set<MolAtom*> > >::const_iterator
 
 2707              chain=pos->mvRotatedChainList.begin();
 
 2708              chain!=pos->mvRotatedChainList.end();++chain)
 
 2710             cout<<
"    -"<<chain->first->GetName()<<
":";
 
 2711             for(set<MolAtom*>::const_iterator pos1=chain->second.begin();
 
 2712                 pos1!=chain->second.end();++pos1)
 
 2713                cout<<(*pos1)->GetName()<<
"  ";
 
 2725          cout<<
"      FLIP ACCEPTED ("<<float(ctflip2)/ctflip1*100<<
"% accepted)"<<endl;
 
 2735       TAU_PROFILE_START(timer1);
 
 2738          static const REAL amp=mBaseRotationAmplitude/(REAL)RAND_MAX;
 
 2742                      ((2.*(REAL)rand()-(REAL)RAND_MAX)*amp*mutationAmplitude*mult,
 
 2743                       (REAL)rand(),(REAL)rand(),(REAL)rand());
 
 2745          mClockOrientation.
Click();
 
 2755       TAU_PROFILE_STOP(timer1);
 
 2760             #if 1 // Move as many atoms as possible 
 2774                   if((*at)->GetX()<xmin) xmin=(*at)->GetX();
 
 2775                   if((*at)->GetX()>xmax) xmax=(*at)->GetX();
 
 2776                   if((*at)->GetY()<ymin) ymin=(*at)->GetY();
 
 2777                   if((*at)->GetY()>ymax) ymax=(*at)->GetY();
 
 2778                   if((*at)->GetZ()<zmin) zmin=(*at)->GetZ();
 
 2779                   if((*at)->GetZ()>zmax) zmax=(*at)->GetZ();
 
 2782                REAL dx=(xmax-xmin)/5.,dy=(ymax-ymin)/5.,dz=(zmax-zmin)/5.;
 
 2786                const REAL xc=xmin+rand()/(REAL)RAND_MAX*(xmax-xmin);
 
 2787                const REAL yc=ymin+rand()/(REAL)RAND_MAX*(ymax-ymin);
 
 2788                const REAL zc=zmin+rand()/(REAL)RAND_MAX*(zmax-zmin);
 
 2789                map<MolAtom*,XYZ> v0;
 
 2790                const REAL ax=-4.*log(2.)/(dx*dx);
 
 2791                const REAL ay=-4.*log(2.)/(dy*dy);
 
 2792                const REAL az=-4.*log(2.)/(dz*dz);
 
 2794                   v0[*at]=
XYZ(exp(ax*((*at)->GetX()-xc)*((*at)->GetX()-xc)),
 
 2795                               exp(ay*((*at)->GetY()-yc)*((*at)->GetY()-yc)),
 
 2796                               exp(az*((*at)->GetZ()-zc)*((*at)->GetZ()-zc)));
 
 2799                map<MolAtom*,XYZ> v0;
 
 2802                std::map<MolAtom*,unsigned long> pushedAtoms;
 
 2803                unsigned long idx=rand()%v0.size();
 
 2805                for(
unsigned int i=0;i<idx;i++) at0++;
 
 2806                const REAL xc=(*at0)->GetX();
 
 2807                const REAL yc=(*at0)->GetY();
 
 2808                const REAL zc=(*at0)->GetZ();
 
 2814                   ux=REAL(rand()-RAND_MAX/2);
 
 2815                   uy=REAL(rand()-RAND_MAX/2);
 
 2816                   uz=REAL(rand()-RAND_MAX/2);
 
 2817                   n=sqrt(ux*ux+uy*uy+uz*uz);
 
 2819                ux=ux/n;uy=uy/n;uz=uz/n;
 
 2820                const REAL a=-4.*log(2.)/(2*2);
 
 2822                   for(map<MolAtom*,unsigned long>::iterator at=pushedAtoms.begin() ;at!=pushedAtoms.end();++at)
 
 2823                      v0[at->first]=
XYZ(ux*exp(a*(at->first->GetX()-xc)*(at->first->GetX()-xc)),
 
 2824                                  uy*exp(a*(at->first->GetY()-yc)*(at->first->GetY()-yc)),
 
 2825                                  uz*exp(a*(at->first->GetZ()-zc)*(at->first->GetZ()-zc)));
 
 2827                   for(map<MolAtom*, unsigned long>::iterator at=pushedAtoms.begin() ;at!=pushedAtoms.end();++at)
 
 2828                      v0[at->first]=
XYZ((at->first->GetX()-xc)*ux*exp(a*(at->first->GetX()-xc)*(at->first->GetX()-xc)),
 
 2829                                        (at->first->GetY()-yc)*uy*exp(a*(at->first->GetY()-yc)*(at->first->GetY()-yc)),
 
 2830                                        (at->first->GetZ()-zc)*uz*exp(a*(at->first->GetZ()-zc)*(at->first->GetZ()-zc)));
 
 2835                                       +this->GetBondAngleList().size()
 
 2836                                       +this->GetDihedralAngleList().size());
 
 2837                map<RigidGroup*,std::pair<XYZ,XYZ> > vr;
 
 2839                                              this->GetBondList(),
 
 2840                                              this->GetBondAngleList(),
 
 2841                                              this->GetDihedralAngleList(),
 
 2844             #else // Move atoms belonging to a MD group 
 2849                for(
unsigned int i=0;i<n;++i)++pos;
 
 2850                map<MolAtom*,XYZ> v0;
 
 2851                for(set<MolAtom*>::iterator at=pos->mvpAtom.begin();at!=pos->mvpAtom.end();++at)
 
 2852                   v0[*at]=
XYZ(rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5);
 
 2855                                     +pos->mvpBondAngle.size()
 
 2856                                     +pos->mvpDihedralAngle.size());
 
 2857                map<RigidGroup*,std::pair<XYZ,XYZ> > vr;
 
 2858                float nrjMult=1.0+mutationAmplitude*0.2;
 
 2859                if((rand()%20)==0) nrjMult=4.0;
 
 2863                                              pos->mvpDihedralAngle,
 
 2871             for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
 
 2883                   (*mode)->CalcDeriv();
 
 2885                   for(map<const MolBond*,REAL>::const_iterator pos=(*mode)->mvpBrokenBond.begin();
 
 2886                       pos!=(*mode)->mvpBrokenBond.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2887                   for(map<const MolBondAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenBondAngle.begin();
 
 2888                       pos!=(*mode)->mvpBrokenBondAngle.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2889                   for(map<const MolDihedralAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenDihedralAngle.begin();
 
 2890                       pos!=(*mode)->mvpBrokenDihedralAngle.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2892                   REAL change=(2.*(REAL)rand()-(REAL)RAND_MAX)/(REAL)RAND_MAX;
 
 2895                   if((*mode)->mLLKDeriv>0)
 
 2897                      change -= 0.3*sqrt(llk);
 
 2898                      if(change<-1) change=-1;
 
 2902                      change += 0.3*sqrt(llk);
 
 2903                      if(change>1) change=1;
 
 2905                   (*mode)->Print(cout);
 
 2906                   change *= mutationAmplitude * (*mode)->mBaseAmplitude;
 
 2907                   cout <<
"      Change="<<change<<
" (dLLK= "<<(*mode)->mLLKDeriv<<
"), llk= "<<llk<<
"     ?->"<<llk+(*mode)->mLLKDeriv*change<<endl;
 
 2909                   (*mode)->Stretch(change);
 
 2912                   for(map<const MolBond*,REAL>::const_iterator pos=(*mode)->mvpBrokenBond.begin();
 
 2913                       pos!=(*mode)->mvpBrokenBond.end();++pos)
 
 2915                       cout<<
"      "<<pos->first->GetName()<<
", llk= "<<pos->first->GetLogLikelihood(
false,
false)
 
 2916                           <<
"   ?->"<<pos->first->GetLogLikelihood(
false,
false)+pos->first->GetDeriv((*mode)->mDerivXYZ,
true)*change
 
 2917                           <<
"?   (deriv="<<pos->first->GetDeriv((*mode)->mDerivXYZ)<<
", "<<pos->first->GetDeriv((*mode)->mDerivXYZ,
true);
 
 2918                       cout<<
")  ->" <<pos->first->GetLogLikelihood()<<endl;
 
 2919                      llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2921                   for(map<const MolBondAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenBondAngle.begin();
 
 2922                       pos!=(*mode)->mvpBrokenBondAngle.end();++pos)
 
 2924                       cout<<
"      "<<pos->first->GetName()<<
", llk= "<<pos->first->GetLogLikelihood(
false,
false)
 
 2925                           <<
"   ?->"<<pos->first->GetLogLikelihood(
false,
false)+pos->first->GetDeriv((*mode)->mDerivXYZ,
true)*change
 
 2926                           <<
"?   (deriv="<<pos->first->GetDeriv((*mode)->mDerivXYZ)<<
", "<<pos->first->GetDeriv((*mode)->mDerivXYZ,
true);
 
 2927                       cout<<
")  ->" <<pos->first->GetLogLikelihood()<<endl;
 
 2928                      llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2930                   for(map<const MolDihedralAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenDihedralAngle.begin();
 
 2931                       pos!=(*mode)->mvpBrokenDihedralAngle.end();++pos)
 
 2933                       cout<<
"      "<<pos->first->GetName()<<
", llk= "<<pos->first->GetLogLikelihood(
false,
false)
 
 2934                           <<
"   ?->"<<pos->first->GetLogLikelihood(
false,
false)+pos->first->GetDeriv((*mode)->mDerivXYZ,
true)*change
 
 2935                           <<
"?   (deriv="<<pos->first->GetDeriv((*mode)->mDerivXYZ)<<
", "<<pos->first->GetDeriv((*mode)->mDerivXYZ,
true);
 
 2936                       cout<<
")  ->" <<pos->first->GetLogLikelihood()<<endl;
 
 2937                      llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2939                   cout <<
" -> "<<llk<<endl;
 
 2944             TAU_PROFILE_START(timer2);
 
 2948                if((rand()%2)==0) (*mode)->RandomStretch(mutationAmplitude);
 
 2950             TAU_PROFILE_STOP(timer2);
 
 2956                TAU_PROFILE_START(timer3);
 
 2957                for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
 
 2963                TAU_PROFILE_STOP(timer3);
 
 2965                TAU_PROFILE_START(timer4);
 
 2973                      (*mode)->CalcDeriv();
 
 2975                      for(map<const MolBond*,REAL>::const_iterator pos=(*mode)->mvpBrokenBond.begin();
 
 2976                          pos!=(*mode)->mvpBrokenBond.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2977                      for(map<const MolBondAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenBondAngle.begin();
 
 2978                          pos!=(*mode)->mvpBrokenBondAngle.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2979                      for(map<const MolDihedralAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenDihedralAngle.begin();
 
 2980                          pos!=(*mode)->mvpBrokenDihedralAngle.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
 
 2981                      REAL change=(2.*(REAL)rand()-(REAL)RAND_MAX)/(REAL)RAND_MAX;
 
 2983                      if((*mode)->mLLKDeriv>0)
 
 2986                         if(change<-1) change=-1;
 
 2991                         if(change>1) change=1;
 
 2993                      if(mutationAmplitude<0.5) change *= mutationAmplitude * (*mode)->mBaseAmplitude;
 
 2994                      else change *= 0.5 * (*mode)->mBaseAmplitude;
 
 2995                      (*mode)->Stretch(change);
 
 3003                TAU_PROFILE_STOP(timer4);
 
 3012                   map<MolAtom*,XYZ> v0;
 
 3013                   for(set<MolAtom*>::iterator at=pos->mvpAtom.begin();at!=pos->mvpAtom.end();++at)
 
 3014                      v0[*at]=
XYZ(rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5);
 
 3016                   const REAL nrj0=20*(pos->mvpBond.size()+pos->mvpBondAngle.size()+pos->mvpDihedralAngle.size());
 
 3017                   map<RigidGroup*,std::pair<XYZ,XYZ> > vr;
 
 3019                                                 (
const vector<MolBond*>) (pos->mvpBond),
 
 3020                                                 (
const vector<MolBondAngle*>) (pos->mvpBondAngle),
 
 3021                                                 (
const vector<MolDihedralAngle*>) (pos->mvpDihedralAngle),
 
 3030             mClockLogLikelihood.
Click();
 
 3037       REAL x0=0,y0=0,z0=0;
 
 3038       for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 3047       for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 3055    VFN_DEBUG_EXIT(
"Molecule::GlobalOptRandomMove()",4)
 
 3060    if(  (mClockLogLikelihood>mClockAtomList)
 
 3061       &&(mClockLogLikelihood>mClockBondList)
 
 3062       &&(mClockLogLikelihood>mClockBondAngleList)
 
 3063       &&(mClockLogLikelihood>mClockDihedralAngleList)
 
 3064       &&(mClockLogLikelihood>mClockAtomPosition)
 
 3066    TAU_PROFILE(
"Molecule::GetLogLikelihood()",
"REAL ()",TAU_DEFAULT);
 
 3068    mClockLogLikelihood.
Click();
 
 3081    for(vector<MolBond*>::const_iterator pos=this->GetBondList().begin();pos!=this->GetBondList().end();++pos)
 
 3082       *p++=(*pos)->GetLength();
 
 3083    for(vector<MolBondAngle*>::const_iterator pos=this->GetBondAngleList().begin();pos!=this->GetBondAngleList().end();++pos)
 
 3084       *p++=(*pos)->GetAngle();
 
 3085    for(vector<MolDihedralAngle*>::const_iterator pos=this->GetDihedralAngleList().begin();pos!=this->GetDihedralAngleList().end();++pos)
 
 3086       *p++=(*pos)->GetAngle();
 
 3094    for(vector<MolBond*>::const_iterator pos=this->GetBondList().begin();pos!=this->GetBondList().end();++pos)
 
 3095       *p++=(*pos)->GetLength0();
 
 3096    for(vector<MolBondAngle*>::const_iterator pos=this->GetBondAngleList().begin();pos!=this->GetBondAngleList().end();++pos)
 
 3097       *p++=(*pos)->GetAngle0();
 
 3098    for(vector<MolDihedralAngle*>::const_iterator pos=this->GetDihedralAngleList().begin();pos!=this->GetDihedralAngleList().end();++pos)
 
 3099       *p++=(*pos)->GetAngle0();
 
 3108    for(vector<MolBond*>::const_iterator pos=this->GetBondList().begin();pos!=this->GetBondList().end();++pos)
 
 3109       *p++=1/((*pos)->GetLengthSigma()* (*pos)->GetLengthSigma()+1e-6);
 
 3110    for(vector<MolBondAngle*>::const_iterator pos=this->GetBondAngleList().begin();pos!=this->GetBondAngleList().end();++pos)
 
 3111       *p++=1/((*pos)->GetAngleSigma()* (*pos)->GetAngleSigma()+1e-6);
 
 3112    for(vector<MolDihedralAngle*>::const_iterator pos=this->GetDihedralAngleList().begin();pos!=this->GetDihedralAngleList().end();++pos)
 
 3113       *p++=1/((*pos)->GetAngleSigma()* (*pos)->GetAngleSigma()+1e-6);
 
 3128    cout<<
"Molecule::TagNewBestConfig()"<<endl;
 
 3130       vector<MolBond*>::const_iterator pos;
 
 3133          cout<<
"BondLength="<<(*pos)->GetLength();
 
 3134          (*pos)->XMLOutput(cout);
 
 3138       vector<MolBondAngle*>::const_iterator pos;
 
 3141          cout<<
"BondAngle="<<(*pos)->GetAngle();
 
 3142          (*pos)->XMLOutput(cout);
 
 3146       vector<MolDihedralAngle*>::const_iterator pos;
 
 3149          cout<<
"DihedralAngle="<<(*pos)->GetAngle();
 
 3150          (*pos)->XMLOutput(cout);
 
 3157       for(set<MolAtom*>::iterator at1=pos->mvpAtom.begin();at1!=pos->mvpAtom.end();++at1)
 
 3159          sprintf(buf,
"%5s : ",(*at1)->GetName().c_str());
 
 3161          for(set<MolAtom*>::iterator at2=at1;at2!=pos->mvpAtom.end();++at2)
 
 3163             if(at1==at2) 
continue;
 
 3164             sprintf(buf,
"%5s(%6.3f),",(*at2)->GetName().c_str(),
GetBondLength(**at1,**at2));
 
 3177    VFN_DEBUG_ENTRY(
"Molecule::GetScatteringComponentList()",3)
 
 3179    VFN_DEBUG_EXIT(
"Molecule::GetScatteringComponentList()",3)
 
 3191    VFN_DEBUG_ENTRY(
"Molecule::POVRayDescription()",3)
 
 3192    const REAL xMin=options.
mXmin; 
const REAL xMax=options.mXmax;
 
 3193    const REAL yMin=options.mYmin; 
const REAL yMax=options.mYmax;
 
 3194    const REAL zMin=options.mZmin; 
const REAL zMax=options.mZmax;
 
 3197       VFN_DEBUG_EXIT(
"Molecule::POVRayDescription():No atom to display !",4)
 
 3206    os << 
"// Description of Molecule :" << this->
GetName() <<endl;
 
 3207    vector<CrystMatrix_REAL> vXYZCoords;
 
 3216          vXYZCoords.push_back(this->
GetCrystal().GetSpaceGroup().
 
 3217                         GetAllSymmetrics(x0,y0,z0,
false,
false,
false));
 
 3220    CrystMatrix_int translate(27,3);
 
 3221    translate=  -1,-1,-1,
 
 3250    CrystVector_REAL xSave,ySave,zSave;
 
 3251    const int nbSymmetrics=vXYZCoords[0].rows();
 
 3253    for(
int i=0;i<nbSymmetrics;i++)
 
 3255       VFN_DEBUG_ENTRY(
"Molecule::POVRayDescription():Symmetric#"<<i,3)
 
 3256       for(
unsigned int j=0;j<
mvpAtom.size();j++)
 
 3258          x(j)=vXYZCoords[j](i,0);
 
 3259          y(j)=vXYZCoords[j](i,1);
 
 3260          z(j)=vXYZCoords[j](i,2);
 
 3266          x(0) = fmod((
float) x(0),(
float)1); 
if(x(0)<0) x(0)+=1.;
 
 3267          y(0) = fmod((
float) y(0),(
float)1); 
if(y(0)<0) y(0)+=1.;
 
 3268          z(0) = fmod((
float) z(0),(
float)1); 
if(z(0)<0) z(0)+=1.;
 
 3272          for(
unsigned int j=1;j<
mvpAtom.size();j++)
 
 3282       for(
int j=0;j<translate.rows();j++)
 
 3284          x += translate(j,0);
 
 3285          y += translate(j,1);
 
 3286          z += translate(j,2);
 
 3291          if(  ((x.min()<xMax) && (x.max()>xMin))
 
 3292             &&((y.min()<yMax) && (y.max()>yMin))
 
 3293             &&((z.min()<zMax) && (z.max()>zMin)))
 
 3295             os<<
"  //Symetric#"<<++ct<<endl;
 
 3296             for(
unsigned int k=0;k<
mvpAtom.size();k++)
 
 3298                isinside(k)=((x(k)>=xMin) && (x(k)<=xMax)) && ((y(k)>=yMin) && (y(k)<=yMax)) && ((z(k)>=zMin) && (z(k)<=zMax));
 
 3299                if(isinside(k)) borderdist(k)=0;
 
 3303                   if(xMin>x(k)) borderdist(k)+=(xMin-x(k))*aa*(xMin-x(k))*aa;
 
 3304                   if(yMin>y(k)) borderdist(k)+=(yMin-y(k))*bb*(yMin-y(k))*bb;
 
 3305                   if(zMin>z(k)) borderdist(k)+=(zMin-z(k))*cc*(zMin-z(k))*cc;
 
 3306                   if(xMax<x(k)) borderdist(k)+=(xMax-x(k))*aa*(xMax-x(k))*aa;
 
 3307                   if(yMax<y(k)) borderdist(k)+=(yMax-y(k))*bb*(yMax-y(k))*bb;
 
 3308                   if(zMax<z(k)) borderdist(k)+=(zMax-z(k))*cc*(zMax-z(k))*cc;
 
 3309                   borderdist(k)=sqrt(borderdist(k));
 
 3315                if((
mvpAtom[k]->IsDummy()) || (fout<0.001)) 
continue;
 
 3316                if(options.
mShowHydrogens==
false && (
mvpAtom[k]->GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)) 
continue;
 
 3341                os << 
"    ObjCrystAtom(" 
 3345                   <<
mvpAtom[k]->GetScatteringPower().GetRadius()/3.0<<
"," 
 3346                   <<
"colour_"+
mvpAtom[k]->GetScatteringPower().GetName()<<
"," 
 3350             for(
unsigned int k=0;k<
mvpBond.size();k++)
 
 3352                if(  (
mvpBond[k]->GetAtom1().IsDummy())
 
 3353                   ||(
mvpBond[k]->GetAtom2().IsDummy()) ) 
continue;
 
 3354                if(options.
mShowHydrogens==
false && (  (
mvpBond[k]->GetAtom1().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)
 
 3355                                                     ||(
mvpBond[k]->GetAtom2().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5))) 
continue;
 
 3356                unsigned long n1,n2;
 
 3358                for(n1=0;n1<
mvpAtom.size();n1++)
 
 3360                for(n2=0;n2<
mvpAtom.size();n2++)
 
 3363                if((isinside(n1)==
false) || (isinside(n2)==
false))
 
 3365                if(fout<0.001) 
continue;
 
 3366                REAL x1=x(n1),y1=y(n1),z1=z(n1),
 
 3367                     x2=x(n2),y2=y(n2),z2=z(n2);
 
 3368                REAL dx=x2-x1,dy=y2-y1,dz=z2-z1;
 
 3369                const REAL r=sqrt(abs(dx*dx+dy*dy+dz*dz))+1e-6;
 
 3370                const REAL r1=
mvpAtom[n1]->GetScatteringPower().GetRadius()/3.0;
 
 3371                const REAL r2=
mvpAtom[n2]->GetScatteringPower().GetRadius()/3.0;
 
 3372                x1+=dx/r*r1*sqrt(abs(1-0.1/r1));
 
 3373                y1+=dy/r*r1*sqrt(abs(1-0.1/r1));
 
 3374                z1+=dz/r*r1*sqrt(abs(1-0.1/r1));
 
 3375                x2-=dx/r*r2*sqrt(abs(1-0.1/r2));
 
 3376                y2-=dy/r*r2*sqrt(abs(1-0.1/r2));
 
 3377                z2-=dz/r*r2*sqrt(abs(1-0.1/r2));
 
 3379                os << 
"    ObjCrystBond(" 
 3380                   <<x1<<
","<<y1<<
","<<z1<< 
"," 
 3381                   <<x2<<
","<<y2<<
","<<z2<< 
"," 
 3383                if(
mvpBond[k]->IsFreeTorsion()) os<<
"colour_freebond,";
 
 3384                else os<<
"colour_nonfreebond,";
 
 3385                os<<f<<
","<<fout<<
")"<<endl;
 
 3392       VFN_DEBUG_EXIT(
"Molecule::POVRayDescription():Symmetric#"<<i,3)
 
 3394    VFN_DEBUG_EXIT(
"Molecule::POVRayDescription()",3)
 
 3399                                const REAL xMin,
const REAL xMax,
 
 3400                                const REAL yMin,
const REAL yMax,
 
 3401                                const REAL zMin,
const REAL zMax,
 
 3402                                const bool displayEnantiomer,
 
 3403                                const bool displayNames,
 
 3404                                const bool hideHydrogens)
const 
 3407    VFN_DEBUG_ENTRY(
"Molecule::GLInitDisplayList()",3)
 
 3410       VFN_DEBUG_EXIT(
"Molecule::GLInitDisplayList():No atom to display !",4)
 
 3414    if(
mvpAtom.size()>200) large=
true;
 
 3416    if(displayEnantiomer==
true) en=-1;
 
 3427    const GLfloat colour0[] = {0.0f, 0.0f, 0.0f, 0.0f};
 
 3428    glMaterialfv(GL_FRONT, GL_SPECULAR,  colour0);
 
 3429    glMaterialfv(GL_FRONT, GL_EMISSION,  colour0);
 
 3430    glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
 
 3431    glPolygonMode(GL_FRONT, GL_FILL);
 
 3433    GLUquadricObj* pQuadric = gluNewQuadric();
 
 3435    if(
true==onlyIndependentAtoms)
 
 3439       vector<MolAtom*>::const_iterator pos;
 
 3443          if((*pos)->IsDummy())
continue;
 
 3444          if(hideHydrogens  && ((*pos)->GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)) 
continue;
 
 3445          const float r=(*pos)->GetScatteringPower().GetColourRGB()[0];
 
 3446          const float g=(*pos)->GetScatteringPower().GetColourRGB()[1];
 
 3447          const float b=(*pos)->GetScatteringPower().GetColourRGB()[2];
 
 3448          const float f=(*pos)->GetOccupancy()*this->
GetOccupancy();
 
 3452                GLfloat colourChar [] = {1.0, 1.0, 1.0, 1.0};
 
 3453                GLfloat colourCharRing [] = {1.0, 1.0, 0.8, 1.0};
 
 3454                if((r>0.8)&&(g>0.8)&&(b>0.8))
 
 3456                   colourChar[0] = 0.5;
 
 3457                   colourChar[1] = 0.5;
 
 3458                   colourChar[2] = 0.5;
 
 3460                if((*pos)->IsInRing())
 
 3461                   glMaterialfv(GL_FRONT, GL_EMISSION,  colourCharRing);
 
 3463                   glMaterialfv(GL_FRONT, GL_EMISSION,  colourChar);
 
 3464                glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
 
 3465                glTranslatef((*pos)->X()*en+xc, (*pos)->Y()+yc, (*pos)->Z()+zc);
 
 3466                crystGLPrint((*pos)->GetName());
 
 3470                const GLfloat colourAtom [] = {r, g, b, f};
 
 3471                glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE,   colourAtom);
 
 3472                glTranslatef((*pos)->X()*en+xc, (*pos)->Y()+yc, (*pos)->Z()+zc);
 
 3473                gluSphere(pQuadric,(*pos)->GetScatteringPower().GetRadius()/3.,20,20);
 
 3480       VFN_DEBUG_ENTRY(
"Molecule::GLInitDisplayList():Show all symmetrics",3)
 
 3482       map<const MolAtom*,unsigned long> rix;
 
 3485             for(vector<MolAtom*>::const_iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 3488       vector<CrystMatrix_REAL> vXYZCoords;
 
 3497             vXYZCoords.push_back(this->
GetCrystal().GetSpaceGroup().
 
 3498                            GetAllSymmetrics(x0,y0,z0,
false,
false,
false));
 
 3501       CrystMatrix_int translate(27,3);
 
 3502       translate=  -1,-1,-1,
 
 3531       CrystVector_REAL xSave,ySave,zSave;
 
 3532       const int nbSymmetrics=vXYZCoords[0].rows();
 
 3533       for(
int i=0;i<nbSymmetrics;i++)
 
 3535          VFN_DEBUG_ENTRY(
"Molecule::GLInitDisplayList():Symmetric#"<<i,3)
 
 3536          for(
unsigned int j=0;j<
mvpAtom.size();j++)
 
 3538             x(j)=vXYZCoords[j](i,0);
 
 3539             y(j)=vXYZCoords[j](i,1);
 
 3540             z(j)=vXYZCoords[j](i,2);
 
 3546             x(0) = fmod((
float) x(0),(
float)1); 
if(x(0)<0) x(0)+=1.;
 
 3547             y(0) = fmod((
float) y(0),(
float)1); 
if(y(0)<0) y(0)+=1.;
 
 3548             z(0) = fmod((
float) z(0),(
float)1); 
if(z(0)<0) z(0)+=1.;
 
 3552             for(
unsigned int j=1;j<
mvpAtom.size();j++)
 
 3562          for(
int j=0;j<translate.rows();j++)
 
 3564             x += translate(j,0);
 
 3565             y += translate(j,1);
 
 3566             z += translate(j,2);
 
 3569             if(  ((x.min()<xMax) && (x.max()>xMin))
 
 3570                &&((y.min()<yMax) && (y.max()>yMin))
 
 3571                &&((z.min()<zMax) && (z.max()>zMin)))
 
 3573                for(
unsigned int k=0;k<
mvpAtom.size();k++)
 
 3575                   isinside(k)=((x(k)>=xMin) && (x(k)<=xMax)) && ((y(k)>=yMin) && (y(k)<=yMax)) && ((z(k)>=zMin) && (z(k)<=zMax));
 
 3576                   if(isinside(k)) borderdist(k)=0;
 
 3580                      if(xMin>x(k)) borderdist(k)+=(xMin-x(k))*aa*(xMin-x(k))*aa;
 
 3581                      if(yMin>y(k)) borderdist(k)+=(yMin-y(k))*bb*(yMin-y(k))*bb;
 
 3582                      if(zMin>z(k)) borderdist(k)+=(zMin-z(k))*cc*(zMin-z(k))*cc;
 
 3583                      if(xMax<x(k)) borderdist(k)+=(xMax-x(k))*aa*(xMax-x(k))*aa;
 
 3584                      if(yMax<y(k)) borderdist(k)+=(yMax-y(k))*bb*(yMax-y(k))*bb;
 
 3585                      if(zMax<z(k)) borderdist(k)+=(zMax-z(k))*cc*(zMax-z(k))*cc;
 
 3586                      borderdist(k)=sqrt(borderdist(k));
 
 3592                   sprintf(ch,
"%d %d %d %s %5.2f %5.2f %5.2f d=%5.2f  fout=%5.3f",i,j,k,
mvpAtom[k]->
GetName().c_str(),x(k),y(k),z(k),borderdist(k),fout);
 
 3593                   VFN_DEBUG_MESSAGE(ch,4)
 
 3596                   if(
mvpAtom[k]->IsDummy()) 
continue;
 
 3597                   if(hideHydrogens  && (
mvpAtom[k]->GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)) 
continue;
 
 3598                   if(fout<0.01) 
continue;
 
 3600                      const float r=
mvpAtom[k]->GetScatteringPower().GetColourRGB()[0];
 
 3601                      const float g=
mvpAtom[k]->GetScatteringPower().GetColourRGB()[1];
 
 3602                      const float b=
mvpAtom[k]->GetScatteringPower().GetColourRGB()[2];
 
 3608                            GLfloat colourChar [] = {1.0, 1.0, 1.0, f*fout};
 
 3609                            GLfloat colourCharRing [] = {1.0, 1.0, 0.8, f*fout};
 
 3610                            if((r>0.8)&&(g>0.8)&&(b>0.8))
 
 3612                               colourChar[0] = 0.5;
 
 3613                               colourChar[1] = 0.5;
 
 3614                               colourChar[2] = 0.5;
 
 3617                               glMaterialfv(GL_FRONT, GL_EMISSION,  colourCharRing);
 
 3619                               glMaterialfv(GL_FRONT, GL_EMISSION,  colourChar);
 
 3620                            glRasterPos3f(x(k)*en, y(k), z(k));
 
 3628                            const GLfloat colourAtom [] = {r*fout, g*fout, b*fout, f*fout};
 
 3629                            glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,colourAtom);
 
 3630                            glTranslatef(x(k)*en, y(k), z(k));
 
 3632                               mvpAtom[k]->GetScatteringPower().GetRadius()/3.,20,20);
 
 3636                            const GLfloat colourAtom [] = {r*fout, g*fout, b*fout, f*fout};
 
 3637                            glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,colourAtom);
 
 3638                            glTranslatef(x(k)*en, y(k), z(k));
 
 3639                            gluSphere(pQuadric,.15,10,10);
 
 3644                if(displayNames==
false)
 
 3648                      for(
unsigned int k=0;k<
mvpBond.size();k++)
 
 3650                         if(  (
mvpBond[k]->GetAtom1().IsDummy())
 
 3651                            ||(
mvpBond[k]->GetAtom2().IsDummy()) ) 
continue;
 
 3652                         if(hideHydrogens  && (  (
mvpBond[k]->GetAtom1().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)
 
 3653                                               ||(
mvpBond[k]->GetAtom2().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5))) 
continue;
 
 3654                         const unsigned long n1=rix[&(
mvpBond[k]->GetAtom1())],
 
 3655                                             n2=rix[&(
mvpBond[k]->GetAtom2())];
 
 3657                         if((isinside(n1)==
false) || (isinside(n2)==
false))
 
 3659                         if(fout<0.01) 
continue;
 
 3661                         const float r1=
mvpBond[k]->GetAtom1().GetScatteringPower().GetColourRGB()[0];
 
 3662                         const float g1=
mvpBond[k]->GetAtom1().GetScatteringPower().GetColourRGB()[1];
 
 3663                         const float b1=
mvpBond[k]->GetAtom1().GetScatteringPower().GetColourRGB()[2];
 
 3665                         const float r2=
mvpBond[k]->GetAtom2().GetScatteringPower().GetColourRGB()[0];
 
 3666                         const float g2=
mvpBond[k]->GetAtom2().GetScatteringPower().GetColourRGB()[1];
 
 3667                         const float b2=
mvpBond[k]->GetAtom2().GetScatteringPower().GetColourRGB()[2];
 
 3669                         const GLfloat colourAtom1 [] = {r1, g1, b1, f1*fout};
 
 3670                         const GLfloat colourAtom2 [] = {r2, g2, b2, f2*fout};
 
 3673                            glBegin(GL_LINE_STRIP);
 
 3674                               glMaterialfv(GL_FRONT, GL_SPECULAR,  colourAtom1);
 
 3675                               glMaterialfv(GL_FRONT, GL_EMISSION,  colourAtom1);
 
 3676                               glMaterialfv(GL_FRONT, GL_SHININESS, colourAtom1);
 
 3677                               glVertex3f(x(n1)*en,y(n1),z(n1));
 
 3678                               glVertex3f((x(n1)+x(n2))/2*en,(y(n1)+y(n2))/2,(z(n1)+z(n2))/2);
 
 3679                               glMaterialfv(GL_FRONT, GL_SPECULAR,  colourAtom2);
 
 3680                               glMaterialfv(GL_FRONT, GL_EMISSION,  colourAtom2);
 
 3681                               glMaterialfv(GL_FRONT, GL_SHININESS, colourAtom2);
 
 3682                               glVertex3f(x(n2)*en,y(n2),z(n2));
 
 3686                         const REAL height= sqrt(abs(  (x(n2)-x(n1))*(x(n2)-x(n1))
 
 3687                                                      +(y(n2)-y(n1))*(y(n2)-y(n1))
 
 3688                                                      +(z(n2)-z(n1))*(z(n2)-z(n1))));
 
 3690                            glTranslatef(x(n1)*en, y(n1), z(n1));
 
 3691                            GLUquadricObj *quadobj = gluNewQuadric();
 
 3692                            glRotatef(180,(x(n2)-x(n1))*en,y(n2)-y(n1),z(n2)-z(n1)+height);
 
 3694                            glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colourAtom1);
 
 3695                            gluCylinder(quadobj,.1,.1,height/2,10,1 );
 
 3696                            gluDeleteQuadric(quadobj);
 
 3698                            glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colourAtom2);
 
 3699                            GLUquadricObj *quadobj2 = gluNewQuadric();
 
 3700                            glTranslatef(0, 0, height/2);
 
 3701                            gluCylinder(quadobj2,.1,.1,height/2,10,1 );
 
 3702                            gluDeleteQuadric(quadobj2);
 
 3709                      for(
unsigned int k=0;k<
mvpBond.size();k++)
 
 3711                         if(  (
mvpBond[k]->GetAtom1().IsDummy())
 
 3712                            ||(
mvpBond[k]->GetAtom2().IsDummy()) ) 
continue;
 
 3713                         if(hideHydrogens  && (  (
mvpBond[k]->GetAtom1().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)
 
 3714                                               ||(
mvpBond[k]->GetAtom2().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5))) 
continue;
 
 3715                         const unsigned long n1=rix[&(
mvpBond[k]->GetAtom1())],
 
 3716                                             n2=rix[&(
mvpBond[k]->GetAtom2())];
 
 3718                         if((isinside(n1)==
false) || (isinside(n2)==
false)) fout=exp(-(borderdist(n1)+borderdist(n2))/2);
 
 3719                         if(fout<0.01) 
continue;
 
 3721                         bool isRigidGroup=
false;
 
 3724                            if(  ((*pos)->find(&(
mvpBond[k]->GetAtom1()))!=(*pos)->end())
 
 3725                               &&((*pos)->find(&(
mvpBond[k]->GetAtom2()))!=(*pos)->end()) )
 
 3732                         const GLfloat colour_bondnonfree[]= { 0.2*fout, .2*fout, .2*fout, f*fout };
 
 3733                         const GLfloat colour_bondrigid[]=   { 0.5*fout, .3*fout, .3*fout, f*fout };
 
 3734                         const GLfloat colour_bondfree[]=    { 0.8*fout, .8*fout, .8*fout, f*fout };
 
 3736                            glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colour_bondrigid);
 
 3739                            if(
mvpBond[k]->IsFreeTorsion())
 
 3740                               glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colour_bondfree);
 
 3742                               glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colour_bondnonfree);
 
 3745                         REAL x1=x(n1),y1=y(n1),z1=z(n1),
 
 3746                         x2=x(n2),y2=y(n2),z2=z(n2);
 
 3747                         REAL dx=x2-x1,dy=y2-y1,dz=z2-z1;
 
 3748                         const REAL r=sqrt(abs(dx*dx+dy*dy+dz*dz))+1e-6;
 
 3749                         const REAL r1=
mvpAtom[n1]->GetScatteringPower().GetRadius()/3.0;
 
 3750                         const REAL r2=
mvpAtom[n2]->GetScatteringPower().GetRadius()/3.0;
 
 3751                         x1+=dx/r*r1*sqrt(abs(1-0.1/r1));
 
 3752                         y1+=dy/r*r1*sqrt(abs(1-0.1/r1));
 
 3753                         z1+=dz/r*r1*sqrt(abs(1-0.1/r1));
 
 3754                         x2-=dx/r*r2*sqrt(abs(1-0.1/r2));
 
 3755                         y2-=dy/r*r2*sqrt(abs(1-0.1/r2));
 
 3756                         z2-=dz/r*r2*sqrt(abs(1-0.1/r2));
 
 3759                            glTranslatef(x1*en, y1, z1);
 
 3760                            GLUquadricObj *quadobj = gluNewQuadric();
 
 3762                            const REAL height= sqrt(abs(  (x2-x1)*(x2-x1)
 
 3765                            glRotatef(180,(x2-x1)*en,y2-y1,z2-z1+height);
 
 3766                            gluCylinder(quadobj,.1,.1,height,10,1 );
 
 3767                            gluDeleteQuadric(quadobj);
 
 3777          VFN_DEBUG_EXIT(
"Molecule::GLInitDisplayList():Symmetric#"<<i,3)
 
 3779       VFN_DEBUG_EXIT(
"Molecule::GLInitDisplayList():Show all symmetrics",3)
 
 3781    gluDeleteQuadric(pQuadric);
 
 3782    VFN_DEBUG_EXIT(
"Molecule::GLInitDisplayList()",3)
 
 3788                        const bool updateDisplay)
 
 3790    VFN_DEBUG_ENTRY(
"Molecule::AddAtom():"<<name,5)
 
 3793    string thename=name;
 
 3794    if(thename==
string(
""))
 
 3797       if(pPow!=0) sprintf(buf,
"%s_%s_%lu",this->
GetName().c_str(),pPow->
GetName().c_str(),
mvpAtom.size()+1);
 
 3798       else sprintf(buf,
"%s_X_%lu",this->
GetName().c_str(),
mvpAtom.size()+1);
 
 3802    mClockAtomPosition.
Click();
 
 3803    mClockAtomScattPow.
Click();
 
 3807                         gpRefParTypeScattConformX,
 
 3808                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 3816                         gpRefParTypeScattConformY,
 
 3817                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 3825                         gpRefParTypeScattConformZ,
 
 3826                         REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 3835    if(molnameasformula || (this->
GetName().size()==0))
 
 3839    VFN_DEBUG_EXIT(
"Molecule::AddAtom()",5)
 
 3844    VFN_DEBUG_ENTRY(
"Molecule::RemoveAtom():"<<atom.GetName(),6)
 
 3845    vector<MolAtom*>::iterator pos=find(
mvpAtom.begin(),
mvpAtom.end(),&atom);
 
 3849                         +
" is not in this Molecule:"+this->
GetName());
 
 3857    for(vector<MolBond*>::iterator posb=
mvpBond.begin();posb!=
mvpBond.end();)
 
 3859       if( (&atom==&((*posb)->GetAtom1())) || (&atom==&((*posb)->GetAtom2())) )
 
 3867       if(  (&atom==&((*posb)->GetAtom1())) || (&atom==&((*posb)->GetAtom2()))
 
 3868          ||(&atom==&((*posb)->GetAtom3())))
 
 3877       if(  (&atom==&((*posb)->GetAtom1())) || (&atom==&((*posb)->GetAtom2()))
 
 3878          ||(&atom==&((*posb)->GetAtom3())) || (&atom==&((*posb)->GetAtom4())))
 
 3882    mClockAtomList.
Click();
 
 3886    if(del) 
delete *pos;
 
 3890    if(molnameasformula || (this->
GetName().size()==0))
 
 3894    VFN_DEBUG_EXIT(
"Molecule::RemoveAtom()",6)
 
 3898 void Molecule::AddNonFlipAtom(
MolAtom &atom)
 
 3900    VFN_DEBUG_ENTRY(
"Molecule::AddNonFlipAtom()",5)
 
 3901    mvNonFlipAtom.push_back(&atom);
 
 3904    mClockFlipGroup.Reset();
 
 3906    VFN_DEBUG_EXIT("
Molecule::AddNonFlipAtom()",5)
 
 3910     for(vector<MolAtom*>::iterator pos=mvNonFlipAtom.begin();pos!=mvNonFlipAtom.end();) {
 
 3911         if(atom.GetName().compare((*pos)->GetName())==0) {
 
 3912             pos = mvNonFlipAtom.erase(pos);
 
 3919 vector<MolAtom*> Molecule::getNonFlipAtomList()
 
 3921     return mvNonFlipAtom;
 
 3924                        const REAL length, 
const REAL sigma, 
const REAL delta,
 
 3925                        const REAL bondOrder,
 
 3926                        const bool updateDisplay)
 
 3928    VFN_DEBUG_ENTRY(
"Molecule::AddBond()",5)
 
 3929    mvpBond.push_back(
new MolBond(atom1,atom2,length,sigma,delta,*
this,bondOrder));
 
 3931    mClockBondList.
Click();
 
 3933    VFN_DEBUG_EXIT(
"Molecule::AddBond()",5)
 
 3938    VFN_DEBUG_ENTRY(
"Molecule::RemoveBond():"<<bond.
GetName(),6)
 
 3939    vector<MolBond*>::iterator pos=find(
mvpBond.begin(),
mvpBond.end(),&bond);
 
 3943                               +
"-"+bond.GetAtom2().GetName()
 
 3944                               +
" is not in this Molecule:"+this->
GetName());
 
 3947    mClockBondList.
Click();
 
 3948    if(del) 
delete *pos;
 
 3951    VFN_DEBUG_EXIT(
"Molecule::RemoveBond():",6)
 
 3957    for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
 
 3959       if(  ((&((*pos)->GetAtom1())==&at1)&&(&((*pos)->GetAtom2())==&at2))
 
 3960          ||((&((*pos)->GetAtom1())==&at2)&&(&((*pos)->GetAtom2())==&at1)))
 
 3967    for(vector<MolBond*>::iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
 
 3969       if(  ((&((*pos)->GetAtom1())==&at1)&&(&((*pos)->GetAtom2())==&at2))
 
 3970          ||((&((*pos)->GetAtom1())==&at2)&&(&((*pos)->GetAtom2())==&at1)))
 
 3977                             const REAL angle, 
const REAL sigma, 
const REAL delta,
 
 3978                             const bool updateDisplay)
 
 3980    VFN_DEBUG_ENTRY(
"Molecule::AddBondAngle()",5)
 
 3983    mClockBondAngleList.
Click();
 
 3985    VFN_DEBUG_EXIT(
"Molecule::AddBondAngle()",5)
 
 3990    VFN_DEBUG_ENTRY(
"Molecule::RemoveBondAngle():"<<angle.GetName(),6)
 
 3994       throw ObjCrystException(
"Molecule::RemoveBondAngle():"+angle.GetAtom1().GetName()
 
 3995                               +
"-"+angle.GetAtom2().GetName()+
"-"+angle.GetAtom3().GetName()
 
 3996                               +
" is not in this Molecule:"+this->
GetName());
 
 3999    mClockBondAngleList.
Click();
 
 4000    if(del) 
delete *pos;
 
 4003    VFN_DEBUG_EXIT(
"Molecule::RemoveBondAngle():",6)
 
 4011    for(vector<MolBondAngle*>::const_iterator pos=
mvpBondAngle.begin();
 
 4014       if(  (&((*pos)->GetAtom2())==&at2)
 
 4015          &&(  ((&((*pos)->GetAtom1())==&at1)&&(&((*pos)->GetAtom3())==&at3))
 
 4016             ||((&((*pos)->GetAtom1())==&at3)&&(&((*pos)->GetAtom3())==&at1))))
 
 4024                                 const REAL angle, 
const REAL sigma, 
const REAL delta,
 
 4025                                 const bool updateDisplay)
 
 4027    VFN_DEBUG_ENTRY(
"Molecule::AddDihedralAngle()",5)
 
 4029                                                    angle,sigma,delta,*
this));
 
 4031    mClockDihedralAngleList.
Click();
 
 4033    VFN_DEBUG_EXIT(
"Molecule::AddDihedralAngle()",5)
 
 4038    VFN_DEBUG_ENTRY(
"Molecule::RemoveDihedralAngle():"<<angle.GetName(),6)
 
 4043       throw ObjCrystException(
"Molecule::RemoveDihedralAngle():"+angle.GetAtom1().GetName()
 
 4044                               +
"-"+angle.GetAtom2().GetName()+
"-"+angle.GetAtom3().GetName()
 
 4045                               +
"-"+angle.GetAtom4().GetName()
 
 4046                               +
" is not in this Molecule:"+this->
GetName());
 
 4049    mClockDihedralAngleList.
Click();
 
 4050    if(del) 
delete *pos;
 
 4053    VFN_DEBUG_ENTRY(
"Molecule::RemoveDihedralAngle():",6)
 
 4061    for(vector<MolDihedralAngle*>::const_iterator pos=
mvpDihedralAngle.begin();
 
 4064       if(  (  ((&((*pos)->GetAtom1())==&at1)&&(&((*pos)->GetAtom2())==&at2))
 
 4065             &&((&((*pos)->GetAtom3())==&at3)&&(&((*pos)->GetAtom4())==&at4)))
 
 4066          ||(  ((&((*pos)->GetAtom4())==&at1)&&(&((*pos)->GetAtom3())==&at2))
 
 4067             &&((&((*pos)->GetAtom2())==&at3)&&(&((*pos)->GetAtom1())==&at4))))
 
 4074                              const bool updateDisplay)
 
 4077   #ifdef RIGID_BODY_STRICT_EXPERIMENTAL 
 4089     sprintf(buf,
"RigidGroup%d_x",i);
 
 4091                     gpRefParTypeScattConformX,
 
 4092                     REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 4099     sprintf(buf,
"RigidGroup%d_y",i);
 
 4101                     gpRefParTypeScattConformY,
 
 4102                     REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 4109     sprintf(buf,
"RigidGroup%d_z",i);
 
 4111                     gpRefParTypeScattConformZ,
 
 4112                     REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
 
 4119     sprintf(buf,
"RigidGroup%d_Q1",i);
 
 4121                     gpRefParTypeScattConform,
 
 4122                     REFPAR_DERIV_STEP_ABSOLUTE,
true,
false,
true,
false,1.,1.);
 
 4129     sprintf(buf,
"RigidGroup%d_Q2",i);
 
 4131                     gpRefParTypeScattConform,
 
 4132                     REFPAR_DERIV_STEP_ABSOLUTE,
true,
false,
true,
false,1.,1.);
 
 4139     sprintf(buf,
"RigidGroup%d_Q3",i);
 
 4141                     gpRefParTypeScattConform,
 
 4142                     REFPAR_DERIV_STEP_ABSOLUTE,
true,
false,
true,
false,1.,1.);
 
 4149   mClockRigidGroup.
Click();
 
 4157    #ifdef RIGID_BODY_STRICT_EXPERIMENTAL 
 4170    if(del) 
delete *pos;
 
 4178 const MolAtom &Molecule::GetAtom(
unsigned int i)
const{
return *
mvpAtom[i];}
 
 4180 MolAtom &Molecule::GetAtom(
const string &name){
return **(this->
FindAtom(name));}
 
 4182 const MolAtom &Molecule::GetAtom(
const string &name)
const{
return **(this->
FindAtom(name));}
 
 4186    VFN_DEBUG_ENTRY(
"Molecule::OptimizeConformation()",5)
 
 4190                                             ANNEALING_EXPONENTIAL,10,.1);
 
 4193    mIsSelfOptimizing=
true;
 
 4194    globalOptObj.
Optimize(nb,
false,stopCost);
 
 4195    mIsSelfOptimizing=
false;
 
 4197    mClockFlipGroup.
Reset();
 
 4198    mClockRotorGroup.
Reset();
 
 4199    VFN_DEBUG_EXIT(
"Molecule::OptimizeConformation()",5)
 
 4205    for(
unsigned i = 0; i < nbStep; ++i)
 
 4209       map<MolAtom*,XYZ> grad;
 
 4210       for(vector<MolAtom*>::iterator pos=this->GetAtomList().begin();
 
 4211          pos!=this->GetAtomList().end();++pos) grad[*pos]=
XYZ(0,0,0);
 
 4214       for(vector<MolBond*>::iterator pos=this->GetBondList().begin();
 
 4215          pos!=this->GetBondList().end();++pos)          (*pos)->CalcGradient(grad);
 
 4217       for(vector<MolBondAngle*>::iterator pos=this->GetBondAngleList().begin();
 
 4218          pos!=this->GetBondAngleList().end();++pos)     (*pos)->CalcGradient(grad);
 
 4220       for(vector<MolDihedralAngle*>::iterator pos=this->GetDihedralAngleList().begin();
 
 4221          pos!=this->GetDihedralAngleList().end();++pos) (*pos)->CalcGradient(grad);
 
 4225       for(map<MolAtom*,XYZ>::const_iterator pos=grad.begin();pos!=grad.end();++pos)
 
 4228          sprintf(buf,
"%10s Grad LLK= (%8.3f %8.3f %8.3f)",
 
 4229                pos->first->GetName().c_str(),pos->second.x,pos->second.y,pos->second.z);
 
 4235       for(map<MolAtom*,XYZ>::const_iterator pos=grad.begin();pos!=grad.end();++pos)
 
 4237          if(abs(pos->second.x)>f) f=abs(pos->second.x);
 
 4238          if(abs(pos->second.y)>f) f=abs(pos->second.y);
 
 4239          if(abs(pos->second.z)>f) f=abs(pos->second.z);
 
 4241       if(f>1e-6) f=maxStep/f;
 
 4248          if((*pos)->size()==0) 
continue; 
 
 4249          REAL dx=0,dy=0,dz=0;
 
 4250          for(set<MolAtom *>::const_iterator at=(*pos)->begin();at!=(*pos)->end();++at)
 
 4259          for(set<MolAtom *>::const_iterator at=(*pos)->begin();at!=(*pos)->end();++at)
 
 4267       for(map<MolAtom*,XYZ>::const_iterator pos=grad.begin();pos!=grad.end();++pos)
 
 4269          pos->first->SetX(pos->first->GetX()-pos->second.x*f);
 
 4270          pos->first->SetY(pos->first->GetY()-pos->second.y*f);
 
 4271          pos->first->SetZ(pos->first->GetZ()-pos->second.z*f);
 
 4279                                        const vector<MolBond*> &vb,
const vector<MolBondAngle*> &va,
 
 4280                                        const vector<MolDihedralAngle*> &vd,
 
 4281                                        map<
RigidGroup*,std::pair<XYZ,XYZ> > &vr, REAL nrj0)
 
 4283    const vector<MolBond*> *pvb=&vb;
 
 4284    const vector<MolBondAngle*> *pva=&va;
 
 4285    const vector<MolDihedralAngle*> *pvd=&vd;
 
 4286    map<RigidGroup*,std::pair<XYZ,XYZ> > *pvr=&vr;
 
 4287    if((pvb->size()==0)&&(pva->size()==0)&&(pvd->size()==0))
 
 4289       pvb=&(this->GetBondList());
 
 4290       pva=&(this->GetBondAngleList());
 
 4291       pvd=&(this->GetDihedralAngleList());
 
 4293          (*pvr)[*pos]=make_pair(
XYZ(0,0,0),
XYZ(0,0,0));
 
 4299    REAL e_v,e_k,v_r=1.0;
 
 4300    for(
unsigned i = 0; i < nbStep; ++i)
 
 4303       map<MolAtom*,XYZ> grad;
 
 4304       for(map<MolAtom*,XYZ>::iterator pos=v0.begin();pos!=v0.end();++pos) grad[pos->first]=
XYZ(0,0,0);
 
 4308       for(vector<MolBond*>::const_iterator pos=pvb->begin();pos!=pvb->end();++pos)
 
 4310          (*pos)->CalcGradient(grad);
 
 4311          e_v+=(*pos)->GetLogLikelihood(
false,
false);
 
 4314       for(vector<MolBondAngle*>::const_iterator pos=pva->begin();pos!=pva->end();++pos)
 
 4316          (*pos)->CalcGradient(grad);
 
 4317          e_v+=(*pos)->GetLogLikelihood(
false,
false);
 
 4320       for(vector<MolDihedralAngle*>::const_iterator pos=pvd->begin();pos!=pvd->end();++pos)
 
 4322          (*pos)->CalcGradient(grad);
 
 4323          e_v+=(*pos)->GetLogLikelihood(
false,
false);
 
 4328       for(map<MolAtom*,XYZ>::const_iterator pos=v0.begin();pos!=v0.end();++pos)
 
 4329          e_k += 0.5*m*(pos->second.x*pos->second.x + pos->second.y*pos->second.y + pos->second.z*pos->second.z);
 
 4331       if(nrj0==0) nrj0=e_k+e_v;
 
 4335          const REAL de=e_k+e_v-nrj0;
 
 4336          if(de<e_k) v_r=sqrt((e_k-de)/e_k);
 
 4342       sprintf(buf,
"(i) LLK + Ek = %10.3f + %10.3f =%10.3f (nrj0=%10.3f)",e_v,e_k,e_v+e_k,nrj0);
 
 4346       for(map<MolAtom*,XYZ>::iterator pos=v0.begin();pos!=v0.end();++pos)
 
 4348          sprintf(buf,
"%10s xyz= (%8.3f %8.3f %8.3f) v= (%8.3f %8.3f %8.3f) m*a= (%8.3f %8.3f %8.3f)",
 
 4349                pos->first->GetName().c_str(),
 
 4350                pos->first->GetX(),pos->first->GetY(),pos->first->GetZ(),
 
 4351                 v0[*pos].x    ,v0[*pos].y     ,v0[*pos].z,
 
 4352                -grad[*pos].x*im,-grad[*pos].y*im,-grad[*pos].z*im);
 
 4358       for(map<MolAtom*,XYZ>::const_iterator pos=v0.begin();pos!=v0.end();++pos)
 
 4360          const XYZ *pa=&(grad[pos->first]);
 
 4361          pos->first->SetX(pos->first->GetX()+pos->second.x*dt*v_r-0.5*im*dt*dt*pa->x);
 
 4362          pos->first->SetY(pos->first->GetY()+pos->second.y*dt*v_r-0.5*im*dt*dt*pa->y);
 
 4363          pos->first->SetZ(pos->first->GetZ()+pos->second.z*dt*v_r-0.5*im*dt*dt*pa->z);
 
 4366       for(map<MolAtom*,XYZ>::iterator pos=v0.begin();pos!=v0.end();++pos)
 
 4368          const XYZ *pa=&(grad[pos->first]);
 
 4369          pos->second.x = v_r*pos->second.x - pa->x*dt*im;
 
 4370          pos->second.y = v_r*pos->second.y - pa->y*dt*im;
 
 4371          pos->second.z = v_r*pos->second.z - pa->z*dt*im;
 
 4376 const vector<MolAtom*>& Molecule::GetAtomList()
const{
return mvpAtom;}
 
 4377 const vector<MolBond*>& Molecule::GetBondList()
const{
return mvpBond;}
 
 4378 const vector<MolBondAngle*>& Molecule::GetBondAngleList()
const{
return mvpBondAngle;}
 
 4379 const vector<MolDihedralAngle*>& Molecule::GetDihedralAngleList()
const{
return mvpDihedralAngle;}
 
 4381 vector<MolAtom*>& Molecule::GetAtomList(){
return mvpAtom;}
 
 4382 vector<MolBond*>& Molecule::GetBondList(){
return mvpBond;}
 
 4383 vector<MolBondAngle*>& Molecule::GetBondAngleList(){
return mvpBondAngle;}
 
 4384 vector<MolDihedralAngle*>& Molecule::GetDihedralAngleList(){
return mvpDihedralAngle;}
 
 4390 const list<StretchModeBondLength>& Molecule::GetStretchModeBondLengthList()
const{
return mvStretchModeBondLength;}
 
 4391 const list<StretchModeBondAngle>& Molecule::GetStretchModeBondAngleList()
const{
return mvStretchModeBondAngle;}
 
 4392 const list<StretchModeTorsion>& Molecule::GetStretchModeTorsionList()
const{
return mvStretchModeTorsion;}
 
 4398                                const set<MolAtom *> &atoms, 
const REAL angle,
 
 4399                                const bool keepCenter)
 
 4401    const REAL vx=at2.X()-at1.X();
 
 4402    const REAL vy=at2.Y()-at1.Y();
 
 4403    const REAL vz=at2.Z()-at1.Z();
 
 4407                                const set<MolAtom *> &atoms, 
const REAL angle,
 
 4408                                const bool keepCenter)
 
 4410    TAU_PROFILE(
"Molecule::RotateAtomGroup(MolAtom&,vx,vy,vz,...)",
"void (...)",TAU_DEFAULT);
 
 4411    if(atoms.size()==0) 
return;
 
 4412    const REAL x0=at.X();
 
 4413    const REAL y0=at.Y();
 
 4414    const REAL z0=at.Z();
 
 4416    if((fabs(vx)+fabs(vy)+fabs(vz))<1e-6) 
return;
 
 4417    REAL dx=0.,dy=0.,dz=0.;
 
 4418    bool keepc=keepCenter;
 
 4422          ||(this->
GetPar(
mXYZ.data()+2).IsFixed())) keepc=
false;
 
 4424    const REAL ca=cos(angle),sa=sin(angle);
 
 4425    const REAL ca1=1-ca;
 
 4426    const REAL vnorm=1/sqrt(vx*vx+vy*vy+vz*vz);
 
 4427    const REAL ux=vx*vnorm,uy=vy*vnorm,uz=vz*vnorm;
 
 4428    const REAL m00=ca+ux*ux*ca1;
 
 4429    const REAL m01=ux*uy*ca1-uz*sa;
 
 4430    const REAL m02=ux*uz*ca1+uy*sa;
 
 4431    const REAL m10=uy*ux*ca1+uz*sa; 
 
 4432    const REAL m11=ca+uy*uy*ca1;
 
 4433    const REAL m12=uy*uz*ca1-ux*sa;
 
 4434    const REAL m20=uz*ux*ca1-uy*sa;
 
 4435    const REAL m21=uz*uy*ca1+ux*sa;
 
 4436    const REAL m22=ca+uz*uz*ca1;
 
 4437    for(set<MolAtom *>::const_iterator pos=atoms.begin();pos!=atoms.end();++pos)
 
 4445       const REAL x=(*pos)->X() - x0;
 
 4446       const REAL y=(*pos)->Y() - y0;
 
 4447       const REAL z=(*pos)->Z() - z0;
 
 4449       (*pos)->X() = m00*x+m01*y+m02*z+x0;
 
 4450       (*pos)->Y() = m10*x+m11*y+m12*z+y0;
 
 4451       (*pos)->Z() = m20*x+m21*y+m22*z+z0;
 
 4461    for(set<MolAtom *>::const_iterator pos=atoms.begin();pos!=atoms.end();++pos)
 
 4472       quat.
RotateVector((*pos)->X(),(*pos)->Y(),(*pos)->Z());
 
 4496    mClockAtomPosition.
Click();
 
 4500                                   const REAL dx,
const REAL dy,
const REAL dz,
 
 4501                                   const bool keepCenter)
 
 4503    for(set<MolAtom *>::const_iterator pos=atoms.begin();pos!=atoms.end();++pos)
 
 4509    bool keepc=keepCenter;
 
 4513          ||(this->
GetPar(
mXYZ.data()+2).IsFixed())) keepc=
false;
 
 4516       const REAL r= (REAL)(atoms.size())/(REAL)(this->
GetNbComponent());
 
 4517       REAL dxc=dx*r,dyc=dy*r,dzc=dz*r;
 
 4523    mClockAtomPosition.
Click();
 
 4528    VFN_DEBUG_ENTRY(
"Molecule::RestraintExport()",5)
 
 4529    os<<
"BondName, IdealLength, Length, log(likelihood)"<<endl;
 
 4530    for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
 
 4531         os <<(*pos)->GetName()
 
 4532            <<
", "<<(*pos)->GetLength0()
 
 4533            <<
", "<<(*pos)->GetLength()
 
 4534            <<
", "<<(*pos)->GetLogLikelihood()<<endl;
 
 4535    os<<
"BondAngle, IdealAngle, Angle, log(likelihood)"<<endl;
 
 4536    for(vector<MolBondAngle*>::const_iterator pos=
mvpBondAngle.begin();
 
 4538         os <<(*pos)->GetName()
 
 4539            <<
", "<<(*pos)->Angle0()*180/M_PI
 
 4540            <<
", "<<(*pos)->GetAngle()*180/M_PI
 
 4541            <<
", "<<(*pos)->GetLogLikelihood()<<endl;
 
 4542    os<<
"DihedralAngle, IdealAngle, Angle, log(likelihood)"<<endl;
 
 4543    for(vector<MolDihedralAngle*>::const_iterator pos=
mvpDihedralAngle.begin();
 
 4545         os <<(*pos)->GetName()
 
 4546            <<
", "<<(*pos)->Angle0()*180/M_PI
 
 4547            <<
", "<<(*pos)->GetAngle()*180/M_PI
 
 4548            <<
", "<<(*pos)->GetLogLikelihood()<<endl;
 
 4549    VFN_DEBUG_EXIT(
"Molecule::RestraintExport()",5)
 
 4553    VFN_DEBUG_ENTRY(
"Molecule::RestraintStatus()",5)
 
 4554    for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
 
 4555       cout <<
"Bond "<<(*pos)->GetName()
 
 4556            <<
", IdealLength="<<(*pos)->GetLength0()
 
 4557            <<
", Length="<<(*pos)->GetLength()
 
 4558            <<
", log(likelihood)="<<(*pos)->GetLogLikelihood()<<endl;
 
 4559    for(vector<MolBondAngle*>::const_iterator pos=
mvpBondAngle.begin();
 
 4561       cout <<
"Bond Angle "<<(*pos)->GetName()
 
 4562            <<
", IdealAngle="<<(*pos)->Angle0()*180/M_PI
 
 4563            <<
", Angle="<<(*pos)->GetAngle()*180/M_PI
 
 4564            <<
", log(likelihood)="<<(*pos)->GetLogLikelihood()<<endl;
 
 4565    for(vector<MolDihedralAngle*>::const_iterator pos=
mvpDihedralAngle.begin();
 
 4567       cout <<
"Dihedral Angle "<<(*pos)->GetName()
 
 4568            <<
", IdealAngle="<<(*pos)->Angle0()*180/M_PI
 
 4569            <<
", Angle="<<(*pos)->GetAngle()*180/M_PI
 
 4570            <<
", log(likelihood)="<<(*pos)->GetLogLikelihood()<<endl;
 
 4571    VFN_DEBUG_EXIT(
"Molecule::RestraintStatus()",5)
 
 4592    for(vector<MolBond*>::iterator bond=
mvpBond.begin();bond!=
mvpBond.end();++bond)
 
 4594       MolAtom* at2=&((*bond)->GetAtom1());
 
 4595       MolAtom* at3=&((*bond)->GetAtom2());
 
 4600          if(*c2==at3) 
continue;
 
 4601          if(
GetBondAngle(**c2,*at2,*at3)<(10 *DEG2RAD)) 
continue;
 
 4602          if(
GetBondAngle(**c2,*at2,*at3)>(180*DEG2RAD)) 
continue;
 
 4607             if((*c3==at2)||(*c3==*c2)) 
continue;
 
 4608             if(
GetBondAngle(*at2,*at3,**c3)<(10 *DEG2RAD)) 
continue;
 
 4609             if(
GetBondAngle(*at2,*at3,**c3)>(180*DEG2RAD)) 
continue;
 
 4624 static const REAL svGaussianIntX[51]
 
 4625    ={0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ,  1.1,  1.2,
 
 4626      1.3,  1.4,  1.5,  1.6,  1.7,  1.8,  1.9,  2. ,  2.1,  2.2,  2.3,
 
 4627      2.4,  2.5,  2.6,  2.7,  2.8,  2.9,  3. ,  3.1,  3.2,  3.3,  3.4,
 
 4628      3.5,  3.6,  3.7,  3.8,  3.9,  4. ,  4.1,  4.2,  4.3,  4.4,  4.5,
 
 4629      4.6,  4.7,  4.8,  4.9,  5. };
 
 4633 static const REAL svGaussianIntY[51]
 
 4634    ={0.00000000e+00,   1.00285768e-01,   2.02498389e-01,   3.08782899e-01,
 
 4635      4.21537858e-01,   5.43577677e-01,   6.78339751e-01,   8.30161516e-01,
 
 4636      1.00466359e+00,   1.20929269e+00,   1.45410564e+00,   1.75292008e+00,
 
 4637      2.12502806e+00,   2.59778353e+00,   3.21056320e+00,   4.02091257e+00,
 
 4638      5.11421527e+00,   6.61911896e+00,   8.73249677e+00,   1.17604260e+01,
 
 4639      1.61864569e+01,   2.27870547e+01,   3.28297946e+01,   4.84189023e+01,
 
 4640      7.31071567e+01,   1.12996748e+02,   1.78751757e+02,   2.89337246e+02,
 
 4641      4.79080996e+02,   8.11232960e+02,   1.40443983e+03,   2.48531490e+03,
 
 4642      4.49461494e+03,   8.30539633e+03,   1.56790580e+04,   3.02354014e+04,
 
 4643      5.95525189e+04,   1.19793234e+05,   2.46080284e+05,   5.16182008e+05,
 
 4644      1.10556243e+06,   2.41765307e+06,   5.39775929e+06,   1.23033271e+07,
 
 4645      2.86288375e+07,   6.80050408e+07,   1.64899866e+08,   4.08157783e+08,
 
 4646      1.03122228e+09,   2.65938808e+09,   7.00012860e+09};
 
 4650 static const CubicSpline sInvGaussianInt(svGaussianIntY,svGaussianIntX,51);
 
 4652 REAL FlatGaussianProba(
const REAL x, 
const REAL sigma, 
const REAL delta)
 
 4654    if(abs(x)<delta) 
return 1.0;
 
 4655    const REAL z=(abs(x)-delta)/sigma;
 
 4659 REAL FlatGaussianIntegral(
const REAL x1,
const REAL x2, 
const REAL sigma, 
const REAL delta)
 
 4661    static const REAL SPI2=0.88622692545275794;
 
 4662    if((x1<=-delta)&&(x2<=-delta)) 
return SPI2*sigma*(erf((x2+delta)/sigma)-erf((x1+delta)/sigma));
 
 4663    if((x1<=-delta)&&(x2<= delta)) 
return (x2+delta)-SPI2*sigma*erf((x1+delta)/sigma);
 
 4664    if(x1<=-delta) 
return 2*delta+SPI2*sigma*(erf((x2-delta)/sigma)-erf((x1+delta)/sigma));
 
 4665    if((x1<= delta)&&(x2<= delta)) 
return x2-x1;
 
 4666    if(x1<= delta) 
return SPI2*sigma*erf((x2-delta)/sigma)+(delta-x1);
 
 4667    return SPI2*sigma*(erf((x2-delta)/sigma)-erf((x1-delta)/sigma));
 
 4670 REAL FlatLorentzianProba(
const REAL x, 
const REAL sigma, 
const REAL delta)
 
 4672    if(abs(x)<delta) 
return 1.0;
 
 4673    const REAL z=(abs(x)-delta)/sigma;
 
 4677 REAL FlatLorentzianIntegral(
const REAL x1,
const REAL x2, 
const REAL sigma, 
const REAL delta)
 
 4679    if((x1<=-delta)&&(x2<=-delta)) 
return atan((x2+delta)/sigma)-atan((x1+delta)/sigma);
 
 4680    if((x1<=-delta)&&(x2<= delta)) 
return (x2+delta)-atan((x1+delta)/sigma);
 
 4681    if(x1<=-delta) 
return 2*delta+atan((x2-delta)/sigma)-atan((x1+delta)/sigma);
 
 4682    if((x1<= delta)&&(x2<= delta)) 
return x2-x1;
 
 4683    if(x1<= delta) 
return atan((x2-delta)/sigma)+(delta-x1);
 
 4684    return atan((x2-delta)/sigma)-atan((x1-delta)/sigma);
 
 4698    REAL r=(REAL)rand()/(REAL)RAND_MAX;
 
 4701       REAL x=x0+amplitude*(2*r-1.0);
 
 4702       if(x> delta)x= delta;
 
 4703       if(x<-delta)x=-delta;
 
 4715       Pmin=FlatLorentzianProba((x0+xmin)/2,sigma,delta);
 
 4716       Pmax=FlatLorentzianProba((x0+xmax)/2,sigma,delta);
 
 4721          for(
int i=0;i<5;i++)
 
 4723             const REAL r=Pmax/Pmin;
 
 4724             xmax=x0+r*amplitude;
 
 4725             Pmax=FlatLorentzianProba((x0+xmax)/2,sigma,delta);
 
 4731          for(
int i=0;i<5;i++)
 
 4733             const REAL r=Pmin/Pmax;
 
 4734             xmin=x0-r*amplitude;
 
 4735             Pmin=FlatLorentzianProba((x0+xmin)/2,sigma,delta);
 
 4748       REAL d=(abs(x0)-delta)/sigma;
 
 4750       d=(1-amplitude*d/2)/(1+amplitude*d/2);
 
 4753          xmin=x0-amplitude*2/(1+d);
 
 4754          xmax=x0+amplitude*2*d/(1+d);
 
 4758          xmax=x0+amplitude*2/(1+d);
 
 4759          xmin=x0-amplitude*2*d/(1+d);
 
 4768       REAL ymin=(abs(xmin)-delta)/sigma;
 
 4770       REAL ymax=(abs(xmax)-delta)/sigma;
 
 4772       const REAL y=ymin+(ymax-ymin)*r;
 
 4773       return -tan(y)*sigma-delta;
 
 4780          REAL p0=atan((abs(xmin)-delta)/sigma)*sigma;
 
 4785             REAL ymin=(abs(xmin)-delta)/sigma;
 
 4787             const REAL y=ymin*(REAL)rand()/(REAL)RAND_MAX;
 
 4788             return -delta-tan(y)*sigma;
 
 4792             return -delta+(REAL)rand()/(REAL)RAND_MAX*(xmax+delta);
 
 4797          REAL p0=atan((abs(xmin)-delta)/sigma)*sigma;
 
 4799          REAL p2=atan((xmax-delta)/sigma)*sigma;
 
 4800          const REAL n=p0+p1+p2;
 
 4803             REAL ymin=(abs(xmin)-delta)/sigma;
 
 4805             const REAL y=ymin*(REAL)rand()/(REAL)RAND_MAX;
 
 4806             const REAL x=-delta-tan(y)*sigma;
 
 4811             const REAL x=-delta+(REAL)rand()/(REAL)RAND_MAX*2*delta;
 
 4815          REAL ymax=(xmax-delta)/sigma;
 
 4817          const REAL y=ymax*(REAL)rand()/(REAL)RAND_MAX;
 
 4818          const REAL x=delta+tan(y)*sigma;
 
 4826          return xmin+r*(xmax-xmin);
 
 4829       const REAL p0=delta-xmin;
 
 4830       const REAL p1=atan((xmax-delta)/sigma)*sigma;
 
 4833          return xmin+(REAL)rand()/(REAL)RAND_MAX*(delta-xmin);
 
 4836       REAL ymax=(xmax-delta)/sigma;
 
 4838       const REAL y=ymax*(REAL)rand()/(REAL)RAND_MAX;
 
 4839       return delta+tan(y)*sigma;
 
 4842    REAL ymin=(xmin-delta)/sigma;
 
 4844    REAL ymax=(xmax-delta)/sigma;
 
 4846    const REAL y=ymin+(ymax-ymin)*r;
 
 4847    return tan(y)*sigma+delta;
 
 4850 void TestLorentzianBiasedRandomMove()
 
 4853    REAL x=0,sigma=0.1,delta=0.5,amplitude=0.05;
 
 4855    for(
long i=0;i<400000;i++)
 
 4888                                       const bool respectRestraint)
 
 4893    const REAL l=sqrt(dx*dx+dy*dy+dz*dz+1e-7);
 
 4895    if(l<1e-6) 
return change;
 
 4896    if(respectRestraint && mode.
mpBond!=0)
 
 4898       const REAL d0=l-mode.
mpBond->GetLength0();
 
 4899       const REAL sigma=mode.
mpBond->GetLengthSigma();
 
 4900       const REAL delta=mode.
mpBond->GetLengthDelta();
 
 4901       const REAL max=delta+sigma*5.0;
 
 4904          REAL d1=d0+(REAL)(2*rand()-RAND_MAX)/(REAL)RAND_MAX*amplitude*0.1;
 
 4905          if(d1> delta)d1= delta;
 
 4906          if(d1<-delta)d1=-delta;
 
 4910       if((d0+change)>max) change=max-d0;
 
 4911       else if((d0+change)<(-max)) change=-max-d0;
 
 4915          cout<<
"BOND LENGTH change("<<change<<
"):" 
 4916              <<mode.
mpAtom0->GetName()<<
"-" 
 4918              <<
"(Restraint="<<l-d0<<
"s"<<sigma<<
"d"<<delta<<
"):" 
 4919              <<l<<
"->"<<l+change<<endl;
 
 4923    else change=(2.*(REAL)rand()-(REAL)RAND_MAX)/(REAL)RAND_MAX*amplitude*0.1;
 
 4932                                      const bool respectRestraint)
 
 4941    const REAL vx=dy10*dz12-dz10*dy12;
 
 4942    const REAL vy=dz10*dx12-dx10*dz12;
 
 4943    const REAL vz=dx10*dy12-dy10*dx12;
 
 4946    if((abs(vx)+abs(vy)+abs(vz))<1e-6) 
return change;
 
 4950       const REAL norm= sqrt( (dx10*dx10+dy10*dy10+dz10*dz10)*(dx12*dx12+dy12*dy12+dz12*dz12)+1e-6);
 
 4951       angle0=(dx10*dx12+dy10*dy12+dz10*dz12)/norm;
 
 4952       if(angle0>=1)  angle0=0;
 
 4955          if(angle0<=-1) angle0=M_PI;
 
 4956          else angle0= acos(angle0);
 
 4959       const REAL a0=angle0-mode.
mpBondAngle->GetAngle0();
 
 4960       const REAL sigma=mode.
mpBondAngle->GetAngleSigma();
 
 4961       const REAL delta=mode.
mpBondAngle->GetAngleDelta();
 
 4964          REAL a1=a0+(REAL)(2*rand()-RAND_MAX)/(REAL)RAND_MAX*amplitude*mode.
mBaseAmplitude;
 
 4965          if(a1> delta)a1= delta;
 
 4966          if(a1<-delta)a1=-delta;
 
 4970       if((a0+change)>(delta+sigma*5.0))       change= delta+sigma*5.0-a0;
 
 4971       else if((a0+change)<(-delta-sigma*5.0)) change=-delta-sigma*5.0-a0;
 
 4975          cout<<
"ANGLE change("<<change*RAD2DEG<<
"):" 
 4976              <<mode.
mpAtom0->GetName()<<
"-" 
 4977              <<mode.
mpAtom1->GetName()<<
"-" 
 4979              <<
"(Restraint="<<(angle0-a0)*RAD2DEG<<
"s"<<sigma*RAD2DEG<<
"d"<<delta*RAD2DEG<<
"):" 
 4980              <<angle0*RAD2DEG<<
"->"<<(angle0+change)*RAD2DEG<<endl;
 
 4984    else change=(2.*(REAL)rand()-(REAL)RAND_MAX)/(REAL)RAND_MAX*mode.
mBaseAmplitude*amplitude;
 
 4989                                          const bool respectRestraint)
 
 4995    if((abs(dx)+abs(dy)+abs(dz))<1e-6) 
return change;
 
 5004          REAL a1=a0+(REAL)(2*rand()-RAND_MAX)/(REAL)RAND_MAX*amplitude*mode.
mBaseAmplitude;
 
 5005          if(a1> delta)a1= delta;
 
 5006          if(a1<-delta)a1=-delta;
 
 5010       if((a0+change)>(delta+sigma*5.0))       change= delta+sigma*5.0-a0;
 
 5011       else if((a0+change)<(-delta-sigma*5.0)) change=-delta-sigma*5.0-a0;
 
 5015          cout<<
"TORSION change (" 
 5016              <<mode.
mpAtom1->GetName()<<
"-"<<mode.
mpAtom2->GetName()<<
"):"<<endl
 
 5017              <<
"     initial angle="<<angle0*RAD2DEG
 
 5019              <<(angle0-a0)*RAD2DEG<<
"s"<<sigma*RAD2DEG<<
"d"<<delta*RAD2DEG<<
"):" 
 5020              <<endl<<
"     New angle:"<<(angle0+change)*RAD2DEG<<
", change="<<change*RAD2DEG<<endl;
 
 5024    else change=(REAL)(2.*rand()-RAND_MAX)/(REAL)RAND_MAX*mode.
mBaseAmplitude*amplitude;
 
 5037    mClockAtomPosition.
Click();
 
 5041 void BuildZMatrixRecursive(
long &z,
const long curr,
 
 5042                            const vector<MolAtom*> &vpAtom,
 
 5043                            const map<
MolAtom *, set<MolAtom *> > &connT,
 
 5044                            vector<MolZAtom> &zmatrix,
 
 5045                            const map<const MolAtom*,long> &vIndex,
 
 5046                            vector<long> &vZIndex,
 
 5047                            vector<long> &vrZIndex)
 
 5049    zmatrix[z].mpPow=&(vpAtom[curr]->GetScatteringPower());
 
 5053       map<MolAtom *, set<MolAtom *> >::const_iterator pConn=connT.find(vpAtom[curr]);
 
 5054       const long nc=pConn->second.size();
 
 5055       vector<long> conn(nc);
 
 5056       vector<long> zconn(nc);
 
 5057       vector<long>::iterator pos=conn.begin();
 
 5058       vector<long>::iterator zpos=zconn.begin();
 
 5059       for(set<MolAtom *>::const_iterator pos1=pConn->second.begin();pos1!=pConn->second.end();++pos1)
 
 5061          *pos = vIndex.find(*pos1)->second;
 
 5062          *zpos = vZIndex[*pos];
 
 5063          cout<<(*pos1)->GetName()<<
"("<<*pos<<
","<<*zpos<<
")"<<endl;
 
 5066       sort(conn.begin(),conn.end());
 
 5067       sort(zconn.begin(),zconn.end());
 
 5071       const long b=zconn[nc-1];
 
 5072       zmatrix[z].mBondAtom=b;
 
 5073       zmatrix[z].mBondLength=
GetBondLength(*vpAtom[vrZIndex[b]],*vpAtom[curr]);
 
 5076          const long a=zmatrix[b].mBondAtom;
 
 5077          zmatrix[z].mBondAngleAtom=a;
 
 5078          zmatrix[z].mBondAngle=
GetBondAngle(*vpAtom[vrZIndex[a]],*vpAtom[vrZIndex[b]],*vpAtom[curr]);
 
 5081             const long d=zmatrix[b].mBondAngleAtom;
 
 5082             zmatrix[z].mDihedralAtom=d;
 
 5083             zmatrix[z].mDihedralAngle=fmod(
GetDihedralAngle(*vpAtom[vrZIndex[d]],*vpAtom[vrZIndex[a]],
 
 5084                                                             *vpAtom[vrZIndex[b]],*vpAtom[curr])+2*M_PI,
 
 5089             zmatrix[z].mDihedralAtom=0;
 
 5090             zmatrix[z].mDihedralAngle=0;
 
 5095          zmatrix[z].mBondAngleAtom=0;
 
 5096          zmatrix[z].mBondAngle=0;
 
 5101       zmatrix[z].mBondAtom=0;
 
 5102       zmatrix[z].mBondLength=0;
 
 5106    for(pos=conn.begin();pos!=conn.end();++pos)
 
 5110          if(vZIndex[*pos]==-1)
 
 5111             BuildZMatrixRecursive(z,*pos,vpAtom,connT,zmatrix,vIndex,vZIndex,vrZIndex);
 
 5121    map<const MolAtom*,long> vIndex;
 
 5124       for(vector<MolAtom*>::const_iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 5130       for(
long i=0;i<n;++i)
 
 5138             for(set<MolAtom *>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 5139                if((vIndex[*pos]<i)&&(vIndex[*pos]>b)) b=vIndex[*pos];
 
 5146                const long a= (b==0)?1 : 
mAsZMatrix[b].mBondAtom;
 
 5157                      for(set<MolAtom *>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 5158                         if((vIndex[*pos]<i) && (vIndex[*pos]!=b) && (vIndex[*pos]>d)) d=vIndex[*pos];
 
 5164                      for(set<MolAtom *>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 5165                         if((vIndex[*pos]<i) && (vIndex[*pos]!=a) && (vIndex[*pos]>d)) d=vIndex[*pos];
 
 5171                      for(set<MolAtom *>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 5172                         if(  (vIndex[*pos]<i)  && (vIndex[*pos]!=a)
 
 5173                            &&(vIndex[*pos]!=b) && (vIndex[*pos]>d)) d=vIndex[*pos];
 
 5177                      for(
long j=0;j<i;++j)
 
 5178                         if((j!=a) &&(j!=b) && (j>d)) d=j;
 
 5194       vector<long> vZIndex(n);
 
 5196       vector<long> vrZIndex(n);
 
 5197       for(
long i=0;i<n;++i)
 
 5223                         const map<
MolAtom *, set<MolAtom *> > &connect,
 
 5224                         list<MolAtom *> &atomlist,
 
 5225                         map<set<MolAtom *>,list<MolAtom *> > &ringlist)
 
 5227    list<MolAtom *>::const_iterator f=find(atomlist.begin(),atomlist.end(),currentAtom);
 
 5228    if(f!=atomlist.end())
 
 5231       cout<<currentAtom->GetName()<<
" was already in the list : ring found !"<<endl;
 
 5232       for(list<MolAtom *>::const_iterator atom=atomlist.begin();atom!=atomlist.end();++atom)
 
 5233          cout<<(*atom)->GetName()<<
" ";
 
 5236       set<MolAtom *> ring1;
 
 5237       list<MolAtom *> ring2;
 
 5238       for(list<MolAtom *>::const_iterator pos=f;pos!=atomlist.end();++pos)
 
 5241          ring2.push_back(*pos);
 
 5243       ringlist.insert(make_pair(ring1,ring2));
 
 5247       atomlist.push_back(currentAtom);
 
 5248       map<MolAtom *,set<MolAtom *> >::const_iterator c=connect.find(currentAtom);
 
 5249       set<MolAtom *>::const_iterator pos;
 
 5250       for(pos=c->second.begin();pos!=c->second.end();++pos)
 
 5252          if(*pos==previousAtom) 
continue;
 
 5255       atomlist.pop_back(); 
 
 5262    if(mClockRingList>mClockConnectivityTable) 
return;
 
 5263    VFN_DEBUG_ENTRY(
"Molecule::BuildRingList()",7)
 
 5264    for(vector<MolAtom*>::const_iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 5265       (*pos)->SetIsInRing(
false);
 
 5266    list<MolAtom *> atomlist;
 
 5268    map<set<MolAtom *>,list<MolAtom *> > ringlist;
 
 5269    for(
unsigned long i=0;i<
mvpAtom.size();i++)
 
 5274    for(map<set<MolAtom *>,list<MolAtom *> >::const_iterator pos0=ringlist.begin();pos0!=ringlist.end();pos0++)
 
 5277       std::list<MolAtom*> *pList=&(
mvRing.back().GetAtomList());
 
 5278       #if 1//def __DEBUG__ 
 5279       cout<<
"Found ring:";
 
 5281       for(list<MolAtom *>::const_iterator atom=pos0->second.begin();atom!=pos0->second.end();++atom)
 
 5283          pList->push_back(*atom);
 
 5284          (*atom)->SetIsInRing(
true);
 
 5285          #if 1//def __DEBUG__ 
 5286          cout<<(*atom)->GetName()<<
" ";
 
 5289       #if 1//def __DEBUG__ 
 5294    cout<<
"Rings found :"<<ringlist.size()<<
", "<<
mvRing.size()<<
" unique."<<endl;
 
 5295    mClockRingList.
Click();
 
 5296    VFN_DEBUG_EXIT(
"Molecule::BuildRingList()",7)
 
 5301    if(  (mClockConnectivityTable>mClockBondList)
 
 5302       &&(mClockConnectivityTable>mClockAtomList)) 
return;
 
 5303    VFN_DEBUG_ENTRY(
"Molecule::BuildConnectivityTable()",5)
 
 5304    TAU_PROFILE(
"Molecule::BuildConnectivityTable()",
"void ()",TAU_DEFAULT);
 
 5306    for(
unsigned long i=0;i<
mvpBond.size();++i)
 
 5314       map<MolAtom *,set<MolAtom *> >::const_iterator pos;
 
 5317          cout<<
"Atom "<<pos->first->GetName()<<
" is connected to atoms: ";
 
 5318          set<MolAtom *>::const_iterator pos1;
 
 5319          for(pos1=pos->second.begin();pos1!=pos->second.end();++pos1)
 
 5321             cout<<(*pos1)->GetName()<<
"  ";
 
 5324          if(pos->second.size()>24)
 
 5325             throw ObjCrystException(
"Molecule: one atom ("+pos->first->GetName()+
") has more than 24 bonds !");
 
 5329    mClockConnectivityTable.
Click();
 
 5330    VFN_DEBUG_EXIT(
"Molecule::BuildConnectivityTable()",5)
 
 5334 mpAtom1(&at1),mpAtom2(&at2),mBaseRotationAmplitude(M_PI*0.04)
 
 5339    if(  (mClockRotorGroup>mClockBondList)
 
 5340       &&(mClockRotorGroup>mClockAtomList)
 
 5341       &&(mClockRotorGroup>mClockBondAngleList)
 
 5342       &&(mClockRotorGroup>mClockDihedralAngleList)) 
return;
 
 5343    VFN_DEBUG_ENTRY(
"Molecule::BuildRotorGroup()",5)
 
 5344    TAU_PROFILE(
"Molecule::BuildRotorGroup()",
"void ()",TAU_DEFAULT);
 
 5351    for(
unsigned long i=0;i<
mvpBond.size();++i)
 
 5355               *
const atom2=&(
mvpBond[i]->GetAtom2());
 
 5356       for(
unsigned int j=1;j<=2;++j)
 
 5358          const set<MolAtom*> *pConn;
 
 5367          for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 5369             if((j==1)&&(*pos==atom2)) 
continue;
 
 5370             if((j==2)&&(*pos==atom1)) 
continue;
 
 5401       for(vector<MolAtom*>::const_iterator atom1=this->GetAtomList().begin();
 
 5402           atom1!=this->GetAtomList().end();++atom1)
 
 5405          vector<MolAtom*>::const_iterator atom2=atom1;
 
 5407          for(;atom2!=this->GetAtomList().end();++atom2)
 
 5409             for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 5411                if(*pos==*atom2) 
continue;
 
 5418                set<MolAtom*>::const_iterator check
 
 5439    for(
unsigned int i=1;i<=3;++i)
 
 5441       list<RotorGroup> *pRotorGroup1;
 
 5448       for(list<RotorGroup>::iterator pos1=pRotorGroup1->begin();
 
 5449           pos1!=pRotorGroup1->end();++pos1)
 
 5451          for(
unsigned int j=i;j<=3;++j)
 
 5453             list<RotorGroup> *pRotorGroup2;
 
 5460             for(list<RotorGroup>::iterator pos2=pRotorGroup2->begin();
 
 5461                 pos2!=pRotorGroup2->end();)
 
 5463                if(pos2==pos1) {++pos2;
continue;}
 
 5464                if((  ((pos1->mpAtom1 == pos2->mpAtom1) && (pos1->mpAtom2 == pos2->mpAtom2))
 
 5465                    ||((pos1->mpAtom2 == pos2->mpAtom1) && (pos1->mpAtom1 == pos2->mpAtom2)))
 
 5466                   &&pos1->mvRotatedAtomList.size() == pos2->mvRotatedAtomList.size())
 
 5469                   for(set<MolAtom*>::const_iterator pos=pos1->mvRotatedAtomList.begin();
 
 5470                       pos!=pos1->mvRotatedAtomList.end();++pos)
 
 5472                      set<MolAtom*>::const_iterator tmp=pos2->mvRotatedAtomList.find(*pos);
 
 5473                      if(tmp == pos2->mvRotatedAtomList.end())
 
 5482                      cout<<
"Identical groups:"<<endl;
 
 5484                          <<pos1->mpAtom1->GetName()<<
"-" 
 5485                          <<pos1->mpAtom2->GetName()<<
" : ";
 
 5486                      for(set<MolAtom*>::iterator pos=pos1->mvRotatedAtomList.begin();
 
 5487                          pos!=pos1->mvRotatedAtomList.end();++pos)
 
 5488                         cout<<(*pos)->GetName()<<
"  ";
 
 5491                          <<pos2->mpAtom1->GetName()<<
"-" 
 5492                          <<pos2->mpAtom2->GetName()<<
" : ";
 
 5493                      for(set<MolAtom*>::iterator pos=pos2->mvRotatedAtomList.begin();
 
 5494                          pos!=pos2->mvRotatedAtomList.end();++pos)
 
 5495                         cout<<(*pos)->GetName()<<
"  ";
 
 5498                      pos2=pRotorGroup2->erase(pos2);
 
 5510    for(
unsigned int i=1;i<=3;++i)
 
 5512       list<RotorGroup> *pRotorGroup1;
 
 5519       for(list<RotorGroup>::iterator pos=pRotorGroup1->begin();
 
 5520           pos!=pRotorGroup1->end();)
 
 5523          for(
unsigned int j=0;j<36;++j)
 
 5525             const REAL angle=(REAL)j*M_PI/36.;
 
 5527                                   pos->mvRotatedAtomList,angle);
 
 5535             case 1: cout<<
"Rotation Group around bond :";
break;
 
 5536             case 2: cout<<
"Rotation Group (single chain) around bond :";
break;
 
 5537             case 3: cout<<
"Rotation Group (internal) between :";
break;
 
 5539          cout <<pos->mpAtom1->GetName()<<
"-" 
 5540              <<pos->mpAtom2->GetName()<<
" : ";
 
 5541          for(set<MolAtom *>::iterator pos1=pos->mvRotatedAtomList.begin();
 
 5542              pos1!=pos->mvRotatedAtomList.end();++pos1)
 
 5543             cout<<(*pos1)->GetName()<<
"  ";
 
 5544          cout<<
"   <d(LLK)>="<< llk/36.;
 
 5548             pos = pRotorGroup1->erase(pos);
 
 5559    for(vector<MolBond*>::iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
 
 5560       (*pos)->SetFreeTorsion(
false);
 
 5564       vector<MolBond*>::iterator  bd=this->
FindBond((*pos->mpAtom1),(*pos->mpAtom2));
 
 5565       if(bd!=
mvpBond.end()) (*bd)->SetFreeTorsion(
true);
 
 5568    mClockRotorGroup.
Click();
 
 5569    VFN_DEBUG_EXIT(
"Molecule::BuildRotorGroup()",5)
 
 5575    if(  (mClockStretchModeBondLength>mClockBondList)
 
 5576       &&(mClockStretchModeBondLength>mClockAtomList)
 
 5577       &&(mClockStretchModeBondLength>mClockBondAngleList)
 
 5578       &&(mClockStretchModeBondLength>mClockDihedralAngleList)) 
return;
 
 5580    VFN_DEBUG_ENTRY(
"Molecule::BuildStretchModeBondLength()",7)
 
 5582    TAU_PROFILE(
"Molecule::BuildStretchModeBondLength()",
"void ()",TAU_DEFAULT);
 
 5583    TAU_PROFILE_TIMER(timer1,
"Molecule::BuildStretchModeBondLength 1",
"", TAU_FIELD);
 
 5584    TAU_PROFILE_TIMER(timer2,
"Molecule::BuildStretchModeBondLength 2",
"", TAU_FIELD);
 
 5585    TAU_PROFILE_TIMER(timer3,
"Molecule::BuildStretchModeBondLength 3",
"", TAU_FIELD);
 
 5586    TAU_PROFILE_TIMER(timer4,
"Molecule::BuildStretchModeBondLength 4",
"", TAU_FIELD);
 
 5587    TAU_PROFILE_TIMER(timer5,
"Molecule::BuildStretchModeBondLength 5",
"", TAU_FIELD);
 
 5591    TAU_PROFILE_START(timer1);
 
 5592    for(
unsigned long i=0;i<
mvpBond.size();++i)
 
 5597       for(
unsigned int j=1;j<=2;++j)
 
 5599          const set<MolAtom*> *pConn;
 
 5614          for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 5616             if((j==1)&&(*pos==atom2)) 
continue;
 
 5617             if((j==2)&&(*pos==atom1)) 
continue;
 
 5623          if( ((j==1)&&(ct2>0)) || ((j==2)&&(ct1>0)) )
 
 5642    TAU_PROFILE_STOP(timer1);
 
 5648       TAU_PROFILE_START(timer5);
 
 5650       for(vector<RigidGroup*>::const_iterator group=
mvRigidGroup.begin();
 
 5654          for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
 
 5655             ct += pos->mvTranslatedAtomList.count(*at);
 
 5656          if((ct>0)&&(ct!=(*group)->size()))
 
 5660             cout<<
"       Breaks Rigid Group:";
 
 5661             for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
 
 5662                cout<<(*at)->GetName()<<
" ";
 
 5670       TAU_PROFILE_STOP(timer5);
 
 5674       unsigned long paramSetRandom[5];
 
 5675       for(
unsigned long i=0;i<5;++i)
 
 5677          for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 5679             (*pos)->SetX(100.*rand()/(REAL) RAND_MAX);
 
 5680             (*pos)->SetY(100.*rand()/(REAL) RAND_MAX);
 
 5681             (*pos)->SetZ(100.*rand()/(REAL) RAND_MAX);
 
 5689       TAU_PROFILE_START(timer2);
 
 5691       pos->mvpBrokenBond.clear();
 
 5692       for(vector<MolBond*>::const_iterator r=
mvpBond.begin();r!=
mvpBond.end();++r)
 
 5695          for(set<MolAtom *>::const_iterator at=pos->mvTranslatedAtomList.begin();
 
 5696              at!=pos->mvTranslatedAtomList.end();++at)
 
 5698             if(*at==&((*r)->GetAtom1())) ct++;
 
 5699             if(*at==&((*r)->GetAtom2())) ct++;
 
 5702          if((ct!=0)&&(ct !=2)) pos->mvpBrokenBond.insert(make_pair(*r,0));
 
 5706          int nb=pos->mvpBrokenBond.size();
 
 5707          if(pos->mpBond!=0) nb -= 1;
 
 5708          if(nb>0) keep=
false;
 
 5712       TAU_PROFILE_STOP(timer2);
 
 5718       TAU_PROFILE_START(timer3);
 
 5720       pos->mvpBrokenBondAngle.clear();
 
 5724          for(set<MolAtom *>::const_iterator at=pos->mvTranslatedAtomList.begin();
 
 5725              at!=pos->mvTranslatedAtomList.end();++at)
 
 5727             if(*at==&((*r)->GetAtom1())) ct++;
 
 5728             if(*at==&((*r)->GetAtom2())) ct++;
 
 5729             if(*at==&((*r)->GetAtom3())) ct++;
 
 5732          if((ct==0)||(ct==3)) broken=
false;
 
 5736             for(
unsigned long i=0;i<5;++i)
 
 5739                pos->CalcDeriv(
false);
 
 5740                (*r)->GetLogLikelihood(
true,
true);
 
 5741                d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 5744             if(abs(d)<=0.01) broken=
false;
 
 5746          if(broken) pos->mvpBrokenBondAngle.insert(make_pair(*r,0));
 
 5750          if(pos->mvpBrokenBondAngle.size()>0) keep=
false;
 
 5754       TAU_PROFILE_STOP(timer3);
 
 5760       TAU_PROFILE_START(timer4);
 
 5762       pos->mvpBrokenDihedralAngle.clear();
 
 5766          for(set<MolAtom *>::const_iterator at=pos->mvTranslatedAtomList.begin();
 
 5767              at!=pos->mvTranslatedAtomList.end();++at)
 
 5769             if(*at==&((*r)->GetAtom1())) ct++;
 
 5770             if(*at==&((*r)->GetAtom2())) ct++;
 
 5771             if(*at==&((*r)->GetAtom3())) ct++;
 
 5772             if(*at==&((*r)->GetAtom4())) ct++;
 
 5775          if((ct==0)||(ct==4)) broken=
false;
 
 5779             for(
unsigned long i=0;i<5;++i)
 
 5782                pos->CalcDeriv(
false);
 
 5783                (*r)->GetLogLikelihood(
true,
true);
 
 5784                d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 5787             if(abs(d)<=0.01) broken=
false;
 
 5789          if(broken) pos->mvpBrokenDihedralAngle.insert(make_pair(*r,0));
 
 5793          if(pos->mvpBrokenDihedralAngle.size()>0) keep=
false;
 
 5797       TAU_PROFILE_STOP(timer4);
 
 5800    for(
unsigned long i=0;i<5;++i) this->
ClearParamSet(paramSetRandom[i]);
 
 5801    #if 1//def __DEBUG__ 
 5802    cout<<
"List of Bond Length stretch modes"<<endl;
 
 5806       cout<<
"   Bond:"<<pos->mpAtom0->GetName()<<
"-"<<pos->mpAtom1->GetName()<<
", Translated Atoms:  ";
 
 5807       for(set<MolAtom*>::const_iterator atom=pos->mvTranslatedAtomList.begin();
 
 5808           atom!=pos->mvTranslatedAtomList.end();++atom)
 
 5810          cout<<(*atom)->GetName()<<
",";
 
 5812       if(pos->mpBond!=0) cout<< 
" ; restrained to length="<<pos->mpBond->GetLength0()
 
 5813                              <<
", sigma="<<pos->mpBond->GetLengthSigma()
 
 5814                              <<
", delta="<<pos->mpBond->GetLengthDelta();
 
 5815       if(pos->mvpBrokenBond.size()>0)
 
 5817          cout<<endl<<
"       Broken bonds:";
 
 5818          for(map<const MolBond*,REAL>::const_iterator bond=pos->mvpBrokenBond.begin();
 
 5819              bond!=pos->mvpBrokenBond.end();++bond)
 
 5820             cout<<bond->first->GetName()<<
", ";
 
 5822       if(pos->mvpBrokenBondAngle.size()>0)
 
 5824          cout<<endl<<
"       Broken bond angles:";
 
 5825          for(map<const MolBondAngle*,REAL>::const_iterator angle=pos->mvpBrokenBondAngle.begin();
 
 5826              angle!=pos->mvpBrokenBondAngle.end();++angle)
 
 5827             cout<<angle->first->GetName()<<
", ";
 
 5829       if(pos->mvpBrokenDihedralAngle.size()>0)
 
 5831          cout<<endl<<
"       Broken dihedral angles:";
 
 5832          for(map<const MolDihedralAngle*,REAL>::const_iterator
 
 5833              angle=pos->mvpBrokenDihedralAngle.begin();
 
 5834              angle!=pos->mvpBrokenDihedralAngle.end();++angle)
 
 5835             cout<<angle->first->GetName()<<
", ";
 
 5840    mClockStretchModeBondLength.
Click();
 
 5841    VFN_DEBUG_EXIT(
"Molecule::BuildStretchModeBondLength()",7)
 
 5847    if(  (mClockStretchModeBondAngle>mClockBondList)
 
 5848       &&(mClockStretchModeBondAngle>mClockAtomList)
 
 5849       &&(mClockStretchModeBondAngle>mClockBondAngleList)
 
 5850       &&(mClockStretchModeBondAngle>mClockDihedralAngleList)) 
return;
 
 5852    VFN_DEBUG_ENTRY(
"Molecule::BuildStretchModeBondAngle()",10)
 
 5854    TAU_PROFILE(
"Molecule::BuildStretchModeBondAngle()",
"void ()",TAU_DEFAULT);
 
 5855    TAU_PROFILE_TIMER(timer1,
"Molecule::BuildStretchModeBondAngle 1",
"", TAU_FIELD);
 
 5856    TAU_PROFILE_TIMER(timer2,
"Molecule::BuildStretchModeBondAngle 2",
"", TAU_FIELD);
 
 5857    TAU_PROFILE_TIMER(timer3,
"Molecule::BuildStretchModeBondAngle 3",
"", TAU_FIELD);
 
 5858    TAU_PROFILE_TIMER(timer4,
"Molecule::BuildStretchModeBondAngle 4",
"", TAU_FIELD);
 
 5859    TAU_PROFILE_TIMER(timer5,
"Molecule::BuildStretchModeBondAngle 5",
"", TAU_FIELD);
 
 5863    TAU_PROFILE_START(timer1);
 
 5864    for(
unsigned long i=0;i<
mvpAtom.size();++i)
 
 5868       if(pConn0->size()<2) 
continue;
 
 5869       for(set<MolAtom*>::const_iterator pos1=pConn0->begin();pos1!=pConn0->end();++pos1)
 
 5871          set<MolAtom*>::const_iterator pos2=pos1;
 
 5873          for(;pos2!=pConn0->end();++pos2)
 
 5875             VFN_DEBUG_MESSAGE(
"Molecule::BuildStretchModeBondAngle():"<<i<<
","<<*pos1<<
","<<*pos2,10)
 
 5880                if(&((*pos)->GetAtom2())==
mvpAtom[i])
 
 5882                   if(  ((&((*pos)->GetAtom1())==*pos1)&&(&((*pos)->GetAtom3())==*pos2))
 
 5883                      ||((&((*pos)->GetAtom1())==*pos2)&&(&((*pos)->GetAtom3())==*pos1)))
 
 5890             for(
unsigned int j=1;j<=2;++j)
 
 5892                const set<MolAtom*> *pConn;
 
 5913                for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 5915                   if(*pos==
mvpAtom[i]) 
continue;
 
 5936                      #if 1//def __DEBUG__ 
 5944                      #if 1//def __DEBUG__ 
 5955    TAU_PROFILE_STOP(timer1);
 
 5960       TAU_PROFILE_START(timer5);
 
 5962       for(vector<RigidGroup*>::const_iterator group=
mvRigidGroup.begin();
 
 5966          for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
 
 5967             ct += pos->mvRotatedAtomList.count(*at);
 
 5971             ct += (*group)->count(pos->mpAtom1);
 
 5972             if(ct!=(*group)->size())
 
 5977                cout<<
"       Breaks Rigid Group:";
 
 5978                for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
 
 5979                   cout<<(*at)->GetName()<<
" ";
 
 5988       TAU_PROFILE_STOP(timer5);
 
 5992       unsigned long paramSetRandom[5];
 
 5993       for(
unsigned long i=0;i<5;++i)
 
 5995          for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 5997             (*pos)->SetX(100.*rand()/(REAL) RAND_MAX);
 
 5998             (*pos)->SetY(100.*rand()/(REAL) RAND_MAX);
 
 5999             (*pos)->SetZ(100.*rand()/(REAL) RAND_MAX);
 
 6007       TAU_PROFILE_START(timer2);
 
 6009       pos->mvpBrokenBond.clear();
 
 6010       for(vector<MolBond*>::const_iterator r=
mvpBond.begin();r!=
mvpBond.end();++r)
 
 6013          for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
 
 6014              at!=pos->mvRotatedAtomList.end();++at)
 
 6016             if(*at==&((*r)->GetAtom1())) ct++;
 
 6017             if(*at==&((*r)->GetAtom2())) ct++;
 
 6021          if((ct==0)||(ct==2)) broken=
false;
 
 6025             for(
unsigned long i=0;i<5;++i)
 
 6028                pos->CalcDeriv(
false);
 
 6029                (*r)->GetLogLikelihood(
true,
true);
 
 6030                d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6033             if(abs(d)<=0.01) broken=
false;
 
 6035          if(broken) pos->mvpBrokenBond.insert(make_pair(*r,0));
 
 6039          if(pos->mvpBrokenBond.size()>0) keep=
false;
 
 6043       TAU_PROFILE_STOP(timer2);
 
 6049       TAU_PROFILE_START(timer3);
 
 6051       pos->mvpBrokenBondAngle.clear();
 
 6055          for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
 
 6056              at!=pos->mvRotatedAtomList.end();++at)
 
 6058             if(*at==&((*r)->GetAtom1())) ct++;
 
 6059             if(*at==&((*r)->GetAtom2())) ct++;
 
 6060             if(*at==&((*r)->GetAtom3())) ct++;
 
 6063          if(ct==0) broken=
false;
 
 6064          if(ct==3) broken=
false;
 
 6068             for(
unsigned long i=0;i<5;++i)
 
 6071                pos->CalcDeriv(
false);
 
 6072                (*r)->GetLogLikelihood(
true,
true);
 
 6073                d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6076             if(abs(d)<=0.01) broken=
false;
 
 6078          if(broken) pos->mvpBrokenBondAngle.insert(make_pair(*r,0));
 
 6082          int nb=pos->mvpBrokenBond.size();
 
 6083          if(pos->mpBondAngle!=0) nb -= 1;
 
 6084          if(nb>0) keep=
false;
 
 6088       TAU_PROFILE_STOP(timer3);
 
 6094       TAU_PROFILE_START(timer4);
 
 6096       pos->mvpBrokenDihedralAngle.clear();
 
 6100          for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
 
 6101              at!=pos->mvRotatedAtomList.end();++at)
 
 6103             if(*at==&((*r)->GetAtom1())) ct++;
 
 6104             if(*at==&((*r)->GetAtom2())) ct++;
 
 6105             if(*at==&((*r)->GetAtom3())) ct++;
 
 6106             if(*at==&((*r)->GetAtom4())) ct++;
 
 6109          if(ct==0) broken=
false;
 
 6110          if(ct==4) broken=
false;
 
 6114             for(
unsigned long i=0;i<5;++i)
 
 6117                pos->CalcDeriv(
false);
 
 6118                (*r)->GetLogLikelihood(
true,
true);
 
 6119                d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6122             if(abs(d)<=0.01) broken=
false;
 
 6124          if(broken) pos->mvpBrokenDihedralAngle.insert(make_pair(*r,0));
 
 6128          if(pos->mvpBrokenDihedralAngle.size()>0) keep=
false;
 
 6132       TAU_PROFILE_STOP(timer4);
 
 6135    for(
unsigned long i=0;i<5;++i) this->
ClearParamSet(paramSetRandom[i]);
 
 6137    cout<<
"List of Bond Angle stretch modes"<<endl;
 
 6141       cout<<
"   Angle:"<<pos->mpAtom0->GetName()<<
"-" 
 6142                        <<pos->mpAtom1->GetName()<<
"-" 
 6143                        <<pos->mpAtom2->GetName()<<
", Rotated Atoms:  ";
 
 6144       for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
 
 6145           atom!=pos->mvRotatedAtomList.end();++atom)
 
 6147          cout<<(*atom)->GetName()<<
",";
 
 6149       if(pos->mpBondAngle!=0) cout<< 
" ; restrained to angle="<<pos->mpBondAngle->GetAngle0()*RAD2DEG
 
 6150                                   <<
"�, sigma="<<pos->mpBondAngle->GetAngleSigma()*RAD2DEG
 
 6151                                   <<
"�, delta="<<pos->mpBondAngle->GetAngleDelta()*RAD2DEG<<
"�";
 
 6152       if(pos->mvpBrokenBond.size()>0)
 
 6154          cout<<endl<<
"       Broken bonds:";
 
 6155          for(map<const MolBond*,REAL>::const_iterator bond=pos->mvpBrokenBond.begin();
 
 6156              bond!=pos->mvpBrokenBond.end();++bond)
 
 6157             cout<<bond->first->GetName()<<
", ";
 
 6159       if(pos->mvpBrokenBondAngle.size()>0)
 
 6161          cout<<endl<<
"       Broken bond angles:";
 
 6162          for(map<const MolBondAngle*,REAL>::const_iterator angle=pos->mvpBrokenBondAngle.begin();
 
 6163              angle!=pos->mvpBrokenBondAngle.end();++angle)
 
 6164             cout<<angle->first->GetName()<<
", ";
 
 6166       if(pos->mvpBrokenDihedralAngle.size()>0)
 
 6168          cout<<endl<<
"       Broken dihedral angles:";
 
 6169          for(map<const MolDihedralAngle*,REAL>::const_iterator
 
 6170              angle=pos->mvpBrokenDihedralAngle.begin();
 
 6171              angle!=pos->mvpBrokenDihedralAngle.end();++angle)
 
 6172             cout<<angle->first->GetName()<<
", ";
 
 6177    mClockStretchModeBondAngle.
Click();
 
 6178    VFN_DEBUG_EXIT(
"Molecule::BuildStretchModeBondAngle()",10)
 
 6184    if(  (mClockStretchModeTorsion>mClockBondList)
 
 6185       &&(mClockStretchModeTorsion>mClockAtomList)
 
 6186       &&(mClockStretchModeTorsion>mClockBondAngleList)
 
 6187       &&(mClockStretchModeTorsion>mClockDihedralAngleList)) 
return;
 
 6189    VFN_DEBUG_ENTRY(
"Molecule::BuildStretchModeTorsion()",7)
 
 6190    TAU_PROFILE(
"Molecule::BuildStretchModeTorsion()",
"void ()",TAU_DEFAULT);
 
 6191    TAU_PROFILE_TIMER(timer1,
"Molecule::BuildStretchModeTorsion 1",
"", TAU_FIELD);
 
 6192    TAU_PROFILE_TIMER(timer2,
"Molecule::BuildStretchModeTorsion 2",
"", TAU_FIELD);
 
 6193    TAU_PROFILE_TIMER(timer3,
"Molecule::BuildStretchModeTorsion 3",
"", TAU_FIELD);
 
 6194    TAU_PROFILE_TIMER(timer4,
"Molecule::BuildStretchModeTorsion 4",
"", TAU_FIELD);
 
 6195    TAU_PROFILE_TIMER(timer5,
"Molecule::BuildStretchModeTorsion 5",
"", TAU_FIELD);
 
 6200    for(
unsigned long i=0;i<
mvpBond.size();++i)
 
 6205       for(
unsigned int j=1;j<=2;++j)
 
 6207          const set<MolAtom*> *pConn;
 
 6214          for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 6216             if((*pos==atom2)||(*pos==atom1)) 
continue;
 
 6228             if(  ((&((*dih)->GetAtom2())==atom1) && (&((*dih)->GetAtom3())==atom2))
 
 6229                ||((&((*dih)->GetAtom3())==atom1) && (&((*dih)->GetAtom2())==atom2)))
 
 6231                const unsigned long ct1=
mvStretchModeTorsion.back().mvRotatedAtomList.count(&((*dih)->GetAtom1()));
 
 6232                const unsigned long ct4=
mvStretchModeTorsion.back().mvRotatedAtomList.count(&((*dih)->GetAtom4()));
 
 6260                if( (  ((mode->mpAtom1==atom1)&&(mode->mpAtom2==atom2))
 
 6261                     ||((mode->mpAtom1==atom2)&&(mode->mpAtom2==atom1)))
 
 6283    for(
unsigned long i=0;i<
mvpBond.size();++i)
 
 6288       for(
unsigned int j=1;j<=2;++j)
 
 6290          const set<MolAtom*> *pConn;
 
 6294          for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 6296             if((*pos==atom2)||(*pos==atom1)) 
continue;
 
 6310                if(  ((&((*dih)->GetAtom2())==atom1) && (&((*dih)->GetAtom3())==atom2))
 
 6311                   ||((&((*dih)->GetAtom3())==atom1) && (&((*dih)->GetAtom2())==atom2)))
 
 6313                   const unsigned long ct1=
mvStretchModeTorsion.back().mvRotatedAtomList.count(&((*dih)->GetAtom1()));
 
 6314                   const unsigned long ct4=
mvStretchModeTorsion.back().mvRotatedAtomList.count(&((*dih)->GetAtom4()));
 
 6342                   if( (  ((mode->mpAtom1==atom1)&&(mode->mpAtom2==atom2))
 
 6343                        ||((mode->mpAtom1==atom2)&&(mode->mpAtom2==atom1)))
 
 6371       TAU_PROFILE_START(timer5);
 
 6373       for(vector<RigidGroup*>::const_iterator group=
mvRigidGroup.begin();
 
 6377          for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
 
 6378             ct += pos->mvRotatedAtomList.count(*at);
 
 6382             ct += (*group)->count(pos->mpAtom1);
 
 6383             ct += (*group)->count(pos->mpAtom2);
 
 6384             if(ct!=(*group)->size())
 
 6389                cout<<
"       Breaks Rigid Group:";
 
 6390                for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
 
 6391                   cout<<(*at)->GetName()<<
" ";
 
 6400       TAU_PROFILE_STOP(timer5);
 
 6404       unsigned long paramSetRandom[5];
 
 6405       for(
unsigned long i=0;i<5;++i)
 
 6407          for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 6409             (*pos)->SetX(100.*rand()/(REAL) RAND_MAX);
 
 6410             (*pos)->SetY(100.*rand()/(REAL) RAND_MAX);
 
 6411             (*pos)->SetZ(100.*rand()/(REAL) RAND_MAX);
 
 6419       TAU_PROFILE_START(timer2);
 
 6421       pos->mvpBrokenBond.clear();
 
 6422       for(vector<MolBond*>::const_iterator r=
mvpBond.begin();r!=
mvpBond.end();++r)
 
 6425          for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
 
 6426              at!=pos->mvRotatedAtomList.end();++at)
 
 6428             if(*at==&((*r)->GetAtom1())) ct++;
 
 6429             if(*at==&((*r)->GetAtom2())) ct++;
 
 6433          if((ct==0)||(ct==2)) broken=
false;
 
 6437             for(
unsigned long i=0;i<5;++i)
 
 6440                pos->CalcDeriv(
false);
 
 6441                (*r)->GetLogLikelihood(
true,
true);
 
 6442                d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6445             if(abs(d)<=0.01) broken=
false;
 
 6447          if(broken) pos->mvpBrokenBond.insert(make_pair(*r,0));
 
 6451          if(pos->mvpBrokenBond.size()>0) keep=
false;
 
 6455       TAU_PROFILE_STOP(timer2);
 
 6461       TAU_PROFILE_START(timer3);
 
 6463       pos->mvpBrokenBondAngle.clear();
 
 6467          for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
 
 6468              at!=pos->mvRotatedAtomList.end();++at)
 
 6470             if(*at==&((*r)->GetAtom1())) ct++;
 
 6471             if(*at==&((*r)->GetAtom2())) ct++;
 
 6472             if(*at==&((*r)->GetAtom3())) ct++;
 
 6475          if((ct==0)||(ct==3)) broken=
false;
 
 6479             for(
unsigned long i=0;i<5;++i)
 
 6482                pos->CalcDeriv(
false);
 
 6483                (*r)->GetLogLikelihood(
true,
true);
 
 6484                d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6487             if(abs(d)<=0.01) broken=
false;
 
 6489          if(broken) pos->mvpBrokenBondAngle.insert(make_pair(*r,0));
 
 6493          if(pos->mvpBrokenBond.size()>0) keep=
false;
 
 6497       TAU_PROFILE_STOP(timer3);
 
 6503       TAU_PROFILE_START(timer4);
 
 6505       pos->mvpBrokenDihedralAngle.clear();
 
 6509          for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
 
 6510              at!=pos->mvRotatedAtomList.end();++at)
 
 6512             if(*at==&((*r)->GetAtom1())) ct++;
 
 6513             if(*at==&((*r)->GetAtom2())) ct++;
 
 6514             if(*at==&((*r)->GetAtom3())) ct++;
 
 6515             if(*at==&((*r)->GetAtom4())) ct++;
 
 6518          if((ct==0)||(ct==4)) broken=
false;
 
 6522             for(
unsigned long i=0;i<5;++i)
 
 6525                pos->CalcDeriv(
false);
 
 6526                (*r)->GetLogLikelihood(
true,
true);
 
 6527                d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6530             if(abs(d)<=0.01) broken=
false;
 
 6532          if(broken) pos->mvpBrokenDihedralAngle.insert(make_pair(*r,0));
 
 6536          int nb=pos->mvpBrokenDihedralAngle.size();
 
 6537          if(pos->mpDihedralAngle!=0) nb -= 1;
 
 6538          if(nb>0) keep=
false;
 
 6542       TAU_PROFILE_STOP(timer4);
 
 6545    for(
unsigned long i=0;i<5;++i) this->
ClearParamSet(paramSetRandom[i]);
 
 6553       cout<<
"   Dihedral Angle:" 
 6554           <<pos->mpAtom1->GetName()<<
"-" 
 6555           <<pos->mpAtom2->GetName()<<
", Rotated Atoms:  ";
 
 6556       for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
 
 6557           atom!=pos->mvRotatedAtomList.end();++atom)
 
 6559          cout<<(*atom)->GetName()<<
",";
 
 6561       if(pos->mpDihedralAngle!=0)
 
 6562          cout<<endl<< 
"      ->restrained by dihedral angle "<<pos->mpDihedralAngle->GetName()
 
 6563              <<
"to :"<<pos->mpDihedralAngle->GetAngle0()*RAD2DEG
 
 6564              <<
"�, sigma="<<pos->mpDihedralAngle->GetAngleSigma()*RAD2DEG
 
 6565              <<
"�, delta="<<pos->mpDihedralAngle->GetAngleDelta()*RAD2DEG<<
"�";
 
 6566       if(pos->mvpBrokenBond.size()>0)
 
 6568          cout<<endl<<
"       Broken bonds:";
 
 6569          for(map<const MolBond*,REAL>::const_iterator bond=pos->mvpBrokenBond.begin();
 
 6570              bond!=pos->mvpBrokenBond.end();++bond)
 
 6571             cout<<bond->first->GetName()<<
", ";
 
 6573       if(pos->mvpBrokenBondAngle.size()>0)
 
 6575          cout<<endl<<
"       Broken bond angles:";
 
 6576          for(map<const MolBondAngle*,REAL>::const_iterator angle=pos->mvpBrokenBondAngle.begin();
 
 6577              angle!=pos->mvpBrokenBondAngle.end();++angle)
 
 6578             cout<<angle->first->GetName()<<
", ";
 
 6580       if(pos->mvpBrokenDihedralAngle.size()>0)
 
 6582          cout<<endl<<
"       Broken dihedral angles:";
 
 6583          for(map<const MolDihedralAngle*,REAL>::const_iterator
 
 6584              angle=pos->mvpBrokenDihedralAngle.begin();
 
 6585              angle!=pos->mvpBrokenDihedralAngle.end();++angle)
 
 6586             cout<<angle->first->GetName()<<
", ";
 
 6591    mClockStretchModeTorsion.
Click();
 
 6592    VFN_DEBUG_EXIT(
"Molecule::BuildStretchModeTorsion()",7)
 
 6598    if(  (mClockStretchModeTwist>mClockBondList)
 
 6599       &&(mClockStretchModeTwist>mClockAtomList)
 
 6600       &&(mClockStretchModeTwist>mClockBondAngleList)
 
 6601       &&(mClockStretchModeTwist>mClockDihedralAngleList)) 
return;
 
 6603    VFN_DEBUG_ENTRY(
"Molecule::BuildStretchModeTwist()",7)
 
 6608    for(vector<MolAtom*>::const_iterator atom1=this->GetAtomList().begin();
 
 6609        atom1!=this->GetAtomList().end();++atom1)
 
 6612       vector<MolAtom*>::const_iterator atom2=atom1;
 
 6614       for(;atom2!=this->GetAtomList().end();++atom2)
 
 6616          for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 6618             if(*pos==*atom2) 
continue;
 
 6625             set<MolAtom*>::const_iterator check
 
 6643                   cout<<
"Rejecting StretchModeTwist ";
mvStretchModeTwist.back().Print(cout);cout<<endl;
 
 6658                      if( (  ((mode->mpAtom1==*atom1)&&(mode->mpAtom2==*atom2))
 
 6659                           ||((mode->mpAtom1==*atom2)&&(mode->mpAtom2==*atom1)))
 
 6663                         cout<<
"Duplicate StretchModeTwist ";
mvStretchModeTwist.back().Print(cout);cout<<endl;
 
 6677                   list<StretchModeTorsion>::reverse_iterator mode;
 
 6680                      if( (  ((mode->mpAtom1==*atom1)&&(mode->mpAtom2==*atom2))
 
 6681                           ||((mode->mpAtom1==*atom2)&&(mode->mpAtom2==*atom1)))
 
 6685                         cout<<
"Duplicate StretchModeTwist (with Torsion) ";
mvStretchModeTwist.back().Print(cout);cout<<endl;
 
 6700       unsigned long paramSetRandom[5];
 
 6701       for(
unsigned long i=0;i<5;++i)
 
 6703          for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
 
 6705             (*pos)->SetX(100.*rand()/(REAL) RAND_MAX);
 
 6706             (*pos)->SetY(100.*rand()/(REAL) RAND_MAX);
 
 6707             (*pos)->SetZ(100.*rand()/(REAL) RAND_MAX);
 
 6715       pos->mvpBrokenBond.clear();
 
 6716       pos->mvpBrokenBondAngle.clear();
 
 6717       pos->mvpBrokenDihedralAngle.clear();
 
 6720       cout<<
"   DerivLLK for Twist mode around:" 
 6721           <<pos->mpAtom1->GetName()<<
"-" 
 6722           <<pos->mpAtom2->GetName()<<
": Moving atoms:";
 
 6723       for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
 
 6724           atom!=pos->mvRotatedAtomList.end();++atom)
 
 6725          cout<<(*atom)->GetName()<<
",";
 
 6728       for(vector<MolBond*>::const_iterator r=
mvpBond.begin();r!=
mvpBond.end();++r)
 
 6730          (*r)->GetLogLikelihood(
true,
true);
 
 6732          for(
unsigned long i=0;i<5;++i)
 
 6735             d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6740             cout<<
"       Bond "<<(*r)->GetName()
 
 6741                 <<
": dLength/da="<<d<<endl;
 
 6743             pos->mvpBrokenBond.insert(make_pair(*r,0.0));
 
 6748          (*r)->GetLogLikelihood(
true,
true);
 
 6750          for(
unsigned long i=0;i<5;++i)
 
 6753             d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6758             cout<<
"       BondAngle:"<<(*r)->GetName()<<
": dAngle/da="<<d<<endl;
 
 6760             pos->mvpBrokenBondAngle.insert(make_pair(*r,0.0));
 
 6767          (*r)->GetLogLikelihood(
true,
true);
 
 6769          for(
unsigned long i=0;i<5;++i)
 
 6772             d += abs((*r)->GetDeriv(pos->mDerivXYZ));
 
 6777             cout<<
"       DihedralAngle:"<<(*r)->GetName()<<
": dAngle/da="<<d<<endl;
 
 6779             pos->mvpBrokenDihedralAngle.insert(make_pair(*r,0.0));
 
 6784       for(vector<RigidGroup*>::const_iterator group=
mvRigidGroup.begin();
 
 6788          for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
 
 6789             ct += pos->mvRotatedAtomList.count(*at);
 
 6794             ct += (*group)->count(pos->mpAtom1);
 
 6795             ct += (*group)->count(pos->mpAtom2);
 
 6796             if(ct!=(*group)->size())
 
 6800                cout<<
"       Breaks Rigid Group:"<<ct<<
"!="<<(*group)->size()<<
":";
 
 6801                for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
 
 6802                   cout<<(*at)->GetName()<<
" ";
 
 6811          if(pos->mvpBrokenBond.size()+pos->mvpBrokenBondAngle.size()+pos->mvpBrokenDihedralAngle.size()) keep=
false;
 
 6816    for(
unsigned long i=0;i<5;++i) this->
ClearParamSet(paramSetRandom[i]);
 
 6824       cout<<
"   Twist mode:" 
 6825           <<pos->mpAtom1->GetName()<<
"-" 
 6826           <<pos->mpAtom2->GetName()<<
", Rotated Atoms:  ";
 
 6827       for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
 
 6828           atom!=pos->mvRotatedAtomList.end();++atom)
 
 6830          cout<<(*atom)->GetName()<<
",";
 
 6832       if(pos->mvpBrokenBond.size()>0)
 
 6834          cout<<endl<<
"       Broken bonds:";
 
 6835          for(map<const MolBond*,REAL>::const_iterator bond=pos->mvpBrokenBond.begin();
 
 6836              bond!=pos->mvpBrokenBond.end();++bond)
 
 6837             cout<<bond->first->GetName()<<
", ";
 
 6839       if(pos->mvpBrokenBondAngle.size()>0)
 
 6841          cout<<endl<<
"       Broken bond angles:";
 
 6842          for(map<const MolBondAngle*,REAL>::const_iterator angle=pos->mvpBrokenBondAngle.begin();
 
 6843              angle!=pos->mvpBrokenBondAngle.end();++angle)
 
 6844             cout<<angle->first->GetName()<<
", ";
 
 6846       if(pos->mvpBrokenDihedralAngle.size()>0)
 
 6848          cout<<endl<<
"       Broken dihedral angles:";
 
 6849          for(map<const MolDihedralAngle*,REAL>::const_iterator
 
 6850              angle=pos->mvpBrokenDihedralAngle.begin();
 
 6851              angle!=pos->mvpBrokenDihedralAngle.end();++angle)
 
 6852             cout<<angle->first->GetName()<<
", ";
 
 6857    mClockStretchModeTwist.
Click();
 
 6858    VFN_DEBUG_EXIT(
"Molecule::BuildStretchModeTwist()",7)
 
 6863    VFN_DEBUG_ENTRY(
"Molecule::TuneGlobalOptimRotationAmplitude()",5)
 
 6864    unsigned long initialConfig=this->
CreateParamSet(
"Initial Configuration");
 
 6865    const unsigned int nbTest=100;
 
 6874    REAL displacement=0;
 
 6876    for(
unsigned int j=0;j<nbTest;j++)
 
 6884       REAL xc=0.,yc=0.,zc=0.;
 
 6887          x0[i]=
mvpAtom[i]->GetX(); xc += x0[i];
 
 6888          y0[i]=
mvpAtom[i]->GetY(); yc += y0[i];
 
 6889          z0[i]=
mvpAtom[i]->GetZ(); zc += z0[i];
 
 6899          const REAL dx10=pos->mpAtom0->GetX()-pos->mpAtom1->GetX();
 
 6900          const REAL dy10=pos->mpAtom0->GetY()-pos->mpAtom1->GetY();
 
 6901          const REAL dz10=pos->mpAtom0->GetZ()-pos->mpAtom1->GetZ();
 
 6902          const REAL dx12=pos->mpAtom2->GetX()-pos->mpAtom1->GetX();
 
 6903          const REAL dy12=pos->mpAtom2->GetY()-pos->mpAtom1->GetY();
 
 6904          const REAL dz12=pos->mpAtom2->GetZ()-pos->mpAtom1->GetZ();
 
 6906          const REAL vx=dy10*dz12-dz10*dy12;
 
 6907          const REAL vy=dz10*dx12-dx10*dz12;
 
 6908          const REAL vz=dx10*dy12-dy10*dx12;
 
 6909          this->
RotateAtomGroup(*(pos->mpAtom1),vx,vy,vz,pos->mvRotatedAtomList,0.01,
false);
 
 6915             pos->mBaseAmplitude+=sqrt(abs(dx*dx+dy*dy+dz*dz));
 
 6917          this->
RotateAtomGroup(*(pos->mpAtom1),vx,vy,vz,pos->mvRotatedAtomList,-0.01,
false);
 
 6922          this->
RotateAtomGroup(*(pos->mpAtom1),*(pos->mpAtom2),pos->mvRotatedAtomList,0.01,
false);
 
 6928             pos->mBaseAmplitude+=sqrt(abs(dx*dx+dy*dy+dz*dz));
 
 6930          this->
RotateAtomGroup(*(pos->mpAtom1),*(pos->mpAtom2),pos->mvRotatedAtomList,-0.01,
false);
 
 6933       for(
unsigned int k=0;k<10;++k)
 
 6936                      (mBaseRotationAmplitude,(REAL)rand(),(REAL)rand(),(REAL)rand());
 
 6946             displacement+=sqrt(abs(dx*dx+dy*dy+dz*dz));
 
 6955       pos->mBaseAmplitude=0.1*0.01/pos->mBaseAmplitude;
 
 6956       if(pos->mBaseAmplitude>.17) pos->mBaseAmplitude=.17;
 
 6958       if((pos->mvpBrokenBond.size()+pos->mvpBrokenBondAngle.size()+pos->mvpBrokenDihedralAngle.size())>0)
 
 6960          pos->mBaseAmplitude/=10;
 
 6963       #if 1// def __DEBUG__ 
 6965           <<pos->mpAtom0->GetName()<<
"-" 
 6966           <<pos->mpAtom1->GetName()<<
"-" 
 6967           <<pos->mpAtom2->GetName()<<
":";
 
 6968       for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
 
 6969           atom!=pos->mvRotatedAtomList.end();++atom) cout<<(*atom)->GetName()<<
",";
 
 6970       cout <<
": base rotation="<<pos->mBaseAmplitude*RAD2DEG<<
" free="<<free<<endl;
 
 6978       pos->mBaseAmplitude=0.1*0.01/pos->mBaseAmplitude;
 
 6979       if(pos->mBaseAmplitude>.17) pos->mBaseAmplitude=.17;
 
 6981       if((pos->mvpBrokenBond.size()+pos->mvpBrokenBondAngle.size()+pos->mvpBrokenDihedralAngle.size())>0)
 
 6983          pos->mBaseAmplitude/=10;
 
 6986       #if 1// def __DEBUG__ 
 6988           <<pos->mpAtom1->GetName()<<
"-" 
 6989           <<pos->mpAtom2->GetName()<<
":";
 
 6990       for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
 
 6991           atom!=pos->mvRotatedAtomList.end();++atom) cout<<(*atom)->GetName()<<
",";
 
 6992       cout <<
": base rotation="<<pos->mBaseAmplitude*RAD2DEG<<
" free="<<free<<endl;
 
 7000       pos->mBaseAmplitude=0.1*0.01/pos->mBaseAmplitude;
 
 7001       if(pos->mBaseAmplitude>.17) pos->mBaseAmplitude=.17;
 
 7003       if((pos->mvpBrokenBond.size()+pos->mvpBrokenBondAngle.size()+pos->mvpBrokenDihedralAngle.size())>0)
 
 7005          pos->mBaseAmplitude/=10;
 
 7008       #if 1// def __DEBUG__ 
 7010           <<pos->mpAtom1->GetName()<<
"-" 
 7011           <<pos->mpAtom2->GetName()<<
":";
 
 7012       for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
 
 7013           atom!=pos->mvRotatedAtomList.end();++atom) cout<<(*atom)->GetName()<<
",";
 
 7014       cout <<
": base rotation="<<pos->mBaseAmplitude*RAD2DEG<<
" free="<<free<<endl;
 
 7020       if(displacement>0) mBaseRotationAmplitude*=0.1/displacement;
 
 7021       if(mBaseRotationAmplitude<(0.02*M_PI/20.))
 
 7023          mBaseRotationAmplitude=0.02*M_PI/20.;
 
 7027       if(mBaseRotationAmplitude>(0.02*M_PI*20.))
 
 7029          mBaseRotationAmplitude=0.02*M_PI*20.;
 
 7037    VFN_DEBUG_EXIT(
"Molecule::TuneGlobalOptimRotationAmplitude()",5)
 
 7043    if(  (mClockFlipGroup>mClockConnectivityTable)
 
 7044       &&(mClockFlipGroup>mClockRigidGroup)) 
return;
 
 7045    VFN_DEBUG_ENTRY(
"Molecule::BuildFlipGroup()",5)
 
 7046    TAU_PROFILE(
"Molecule::BuildFlipGroup()",
"void ()",TAU_DEFAULT);
 
 7049    for(vector<MolAtom*>::const_iterator atom0=this->GetAtomList().begin();
 
 7050        atom0!=this->GetAtomList().end();++atom0)
 
 7053       if(pConn->size()<3) 
continue;
 
 7055       for(set<MolAtom*>::const_iterator pos1=pConn->begin();pos1!=pConn->end();++pos1)
 
 7057          for(set<MolAtom*>::const_iterator pos2=pos1;pos2!=pConn->end();++pos2)
 
 7059             if(*pos2==*pos1) 
continue;
 
 7063                bool foundRing=
false;
 
 7064                for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
 
 7066                   if((pos==pos1)||(pos==pos2)) 
continue;
 
 7068                      make_pair(*pos,set<MolAtom*>()));
 
 7069                   mvFlipGroup.back().mvRotatedChainList.back().second.insert(*atom0);
 
 7071                                            mvFlipGroup.back().mvRotatedChainList.back().second);
 
 7072                   mvFlipGroup.back().mvRotatedChainList.back().second.erase(*atom0);
 
 7073                   set<MolAtom*>::const_iterator ringdetect1,ringdetect2;
 
 7074                   ringdetect1=find(
mvFlipGroup.back().mvRotatedChainList.back().second.begin(),
 
 7075                                    mvFlipGroup.back().mvRotatedChainList.back().second.end(),
 
 7077                   ringdetect2=find(
mvFlipGroup.back().mvRotatedChainList.back().second.begin(),
 
 7078                                    mvFlipGroup.back().mvRotatedChainList.back().second.end(),
 
 7080                   if(  (ringdetect1!=
mvFlipGroup.back().mvRotatedChainList.back().second.end())
 
 7081                      ||(ringdetect2!=
mvFlipGroup.back().mvRotatedChainList.back().second.end()))
 
 7084                unsigned long flipSize=0;
 
 7085                for(list<pair<
const MolAtom *,set<MolAtom*> > >::const_iterator
 
 7086                    chain=
mvFlipGroup.back().mvRotatedChainList.begin();
 
 7087                    chain!=
mvFlipGroup.back().mvRotatedChainList.end();++chain)
 
 7088                    flipSize+=chain->second.size();
 
 7095                make_pair(*atom0,set<MolAtom*>()));
 
 7096                mvFlipGroup.back().mvRotatedChainList.back().second.insert(*atom0);
 
 7098                                      mvFlipGroup.back().mvRotatedChainList.back().second);
 
 7100                                      mvFlipGroup.back().mvRotatedChainList.back().second);
 
 7101             mvFlipGroup.back().mvRotatedChainList.back().second.erase(*atom0);
 
 7111       set<MolAtom*> fullset;
 
 7112       for(list<pair<
const MolAtom *,set<MolAtom *> > >::iterator chain=pos->mvRotatedChainList.begin();
 
 7113             chain!=pos->mvRotatedChainList.end();++chain)
 
 7114          for(set<MolAtom *>::const_iterator at=chain->second.begin();at!=chain->second.end();++at)
 
 7115             fullset.insert(*at);
 
 7121          for(set<MolAtom *>::const_iterator at=fullset.begin();at!=fullset.end();++at)
 
 7122             ct+=(*group)->count(*at);
 
 7124          if((ct>0)&&(ct<(*group)->size())) {keep=
false; 
break;}
 
 7129          cout <<
"EXCLUDING flip group (breaking a rigid group): " 
 7130               <<pos->mpAtom0->GetName()<<
",exchanging bonds with " 
 7131               <<pos->mpAtom1->GetName()<<
" and " 
 7132               <<pos->mpAtom2->GetName()<<
", resulting in a 180deg rotation of atoms : ";
 
 7133          for(set<MolAtom*>::iterator pos1=pos->mvRotatedChainList.begin()->second.begin();
 
 7134              pos1!=pos->mvRotatedChainList.begin()->second.end();++pos1)
 
 7135             cout<<(*pos1)->GetName()<<
"  ";
 
 7145        for(
size_t i = 0; i < mvNonFlipAtom.size(); i++) {
 
 7146            if(pos->mpAtom0->GetName().compare(mvNonFlipAtom[i]->
GetName())==0) {
 
 7152            cout <<
"EXCLUDING flip group (central atom is in the non-flip list)"<<endl;
 
 7160    #if 1//def __DEBUG__ 
 7162    for(list<FlipGroup>::iterator pos=
mvFlipGroup.begin();
 
 7165       if(pos->mvRotatedChainList.begin()->first==pos->mpAtom0)
 
 7167          cout <<
"Flip group from atom " 
 7168               <<pos->mpAtom0->GetName()<<
",exchanging bonds with " 
 7169               <<pos->mpAtom1->GetName()<<
" and " 
 7170               <<pos->mpAtom2->GetName()<<
", resulting in a 180� rotation of atoms : ";
 
 7171          for(set<MolAtom*>::iterator pos1=pos->mvRotatedChainList.begin()->second.begin();
 
 7172              pos1!=pos->mvRotatedChainList.begin()->second.end();++pos1)
 
 7173             cout<<(*pos1)->GetName()<<
"  ";
 
 7177          cout <<
"Flip group with respect to: " 
 7178               <<pos->mpAtom1->GetName()<<
"-" 
 7179               <<pos->mpAtom0->GetName()<<
"-" 
 7180               <<pos->mpAtom2->GetName()<<
" : ";
 
 7181          for(list<pair<
const MolAtom *,set<MolAtom*> > >::const_iterator
 
 7182              chain=pos->mvRotatedChainList.begin();
 
 7183              chain!=pos->mvRotatedChainList.end();++chain)
 
 7185             cout<<
"    -"<<chain->first->GetName()<<
":";
 
 7186             for(set<MolAtom*>::const_iterator pos1=chain->second.begin();
 
 7187                 pos1!=chain->second.end();++pos1)
 
 7188                cout<<(*pos1)->GetName()<<
"  ";
 
 7200          cout <<
" -> NOT a free flip, d(llk)="<<dllk;
 
 7203       else cout <<
" -> free flip, d(llk)="<<dllk;
 
 7209    mClockFlipGroup.
Click();
 
 7210    VFN_DEBUG_EXIT(
"Molecule::BuildFlipGroup()",5)
 
 7216    map<const MolBond*,         set<const StretchMode*> > vpBond;
 
 7220       if(mode->mpBond!=0) vpBond[mode->mpBond].insert(&(*mode));
 
 7221       for(map<const MolBond*,REAL>::const_iterator pos=mode->mvpBrokenBond.begin();
 
 7222           pos!=mode->mvpBrokenBond.end();++pos)
 
 7223           vpBond[pos->first].insert(&(*mode));
 
 7227       for(map<const MolBond*,REAL>::const_iterator pos=mode->mvpBrokenBond.begin();
 
 7228           pos!=mode->mvpBrokenBond.end();++pos)
 
 7229           vpBond[pos->first].insert(&(*mode));
 
 7233       for(map<const MolBond*,REAL>::const_iterator pos=mode->mvpBrokenBond.begin();
 
 7234           pos!=mode->mvpBrokenBond.end();++pos)
 
 7235           vpBond[pos->first].insert(&(*mode));
 
 7238       for(map<const MolBond*,REAL>::const_iterator pos=mode->mvpBrokenBond.begin();
 
 7239           pos!=mode->mvpBrokenBond.end();++pos)
 
 7240           vpBond[pos->first].insert(&(*mode));
 
 7244    map<const MolBondAngle*,    set<const StretchMode*> > vpAngle;
 
 7247       for(map<const MolBondAngle*,REAL>::const_iterator pos=mode->mvpBrokenBondAngle.begin();
 
 7248           pos!=mode->mvpBrokenBondAngle.end();++pos)
 
 7249           vpAngle[pos->first].insert(&(*mode));
 
 7254       if(mode->mpBondAngle!=0) vpAngle[mode->mpBondAngle].insert(&(*mode));
 
 7255       for(map<const MolBondAngle*,REAL>::const_iterator pos=mode->mvpBrokenBondAngle.begin();
 
 7256           pos!=mode->mvpBrokenBondAngle.end();++pos)
 
 7257           vpAngle[pos->first].insert(&(*mode));
 
 7261       for(map<const MolBondAngle*,REAL>::const_iterator pos=mode->mvpBrokenBondAngle.begin();
 
 7262           pos!=mode->mvpBrokenBondAngle.end();++pos)
 
 7263           vpAngle[pos->first].insert(&(*mode));
 
 7266       for(map<const MolBondAngle*,REAL>::const_iterator pos=mode->mvpBrokenBondAngle.begin();
 
 7267           pos!=mode->mvpBrokenBondAngle.end();++pos)
 
 7268           vpAngle[pos->first].insert(&(*mode));
 
 7272    map<const MolDihedralAngle*,set<const StretchMode*> > vpDihed;
 
 7275       for(map<const MolDihedralAngle*,REAL>::const_iterator pos=mode->mvpBrokenDihedralAngle.begin();
 
 7276           pos!=mode->mvpBrokenDihedralAngle.end();++pos)
 
 7277           vpDihed[pos->first].insert(&(*mode));
 
 7281       for(map<const MolDihedralAngle*,REAL>::const_iterator pos=mode->mvpBrokenDihedralAngle.begin();
 
 7282           pos!=mode->mvpBrokenDihedralAngle.end();++pos)
 
 7283           vpDihed[pos->first].insert(&(*mode));
 
 7288       if(mode->mpDihedralAngle!=0) vpDihed[mode->mpDihedralAngle].insert(&(*mode));
 
 7289       for(map<const MolDihedralAngle*,REAL>::const_iterator pos=mode->mvpBrokenDihedralAngle.begin();
 
 7290           pos!=mode->mvpBrokenDihedralAngle.end();++pos)
 
 7291           vpDihed[pos->first].insert(&(*mode));
 
 7295       for(map<const MolDihedralAngle*,REAL>::const_iterator pos=mode->mvpBrokenDihedralAngle.begin();
 
 7296           pos!=mode->mvpBrokenDihedralAngle.end();++pos)
 
 7297           vpDihed[pos->first].insert(&(*mode));
 
 7299    for(map<
const MolBond*,set<const StretchMode*> >::const_iterator pos=vpBond.begin();pos!=vpBond.end();++pos)
 
 7301       cout<<
"Bond "<<pos->first->GetName()<<
" is modified by the stretch modes:"<<endl;
 
 7302       for(set<const StretchMode*>::const_iterator mode=pos->second.begin();mode!=pos->second.end();++mode)
 
 7305          (*mode)->Print(cout);
 
 7309    for(map<
const MolBondAngle*,set<const StretchMode*> >::const_iterator pos=vpAngle.begin();pos!=vpAngle.end();++pos)
 
 7311       cout<<
"Bond Angle "<<pos->first->GetName()<<
" is modified by the stretch modes:"<<endl;
 
 7312       for(set<const StretchMode*>::const_iterator mode=pos->second.begin();mode!=pos->second.end();++mode)
 
 7315          (*mode)->Print(cout);
 
 7319    for(map<
const MolDihedralAngle*,set<const StretchMode*> >::const_iterator pos=vpDihed.begin();pos!=vpDihed.end();++pos)
 
 7321       cout<<
"Dihedral Angle "<<pos->first->GetName()<<
" is modified by the stretch modes:"<<endl;
 
 7322       for(set<const StretchMode*>::const_iterator mode=pos->second.begin();mode!=pos->second.end();++mode)
 
 7325          (*mode)->Print(cout);
 
 7336       int nb=mode->mvpBrokenDihedralAngle.size()+mode->mvpBrokenBondAngle.size()+mode->mvpBrokenBond.size();
 
 7337       if(mode->mpBond!=0) nb -= 1;
 
 7345       int nb=mode->mvpBrokenDihedralAngle.size()+mode->mvpBrokenBondAngle.size()+mode->mvpBrokenBond.size();
 
 7346       if(mode->mpBondAngle!=0) nb -= 1;
 
 7354       int nb=mode->mvpBrokenDihedralAngle.size()+mode->mvpBrokenBondAngle.size()+mode->mvpBrokenBond.size();
 
 7355       if(mode->mpDihedralAngle!=0) nb -= 1;
 
 7362       if(mode->mvpBrokenDihedralAngle.size()+mode->mvpBrokenBondAngle.size()+mode->mvpBrokenBond.size()==0)
 
 7372    map<MolAtom*,set<MolAtom*> > vBoundAtoms;
 
 7374    for(vector<MolAtom*>::iterator pos=this->GetAtomList().begin();pos!=this->GetAtomList().end();++pos)
 
 7377    for(vector<MolAtom*>::iterator pat1=this->GetAtomList().begin();pat1!=this->GetAtomList().end();++pat1)
 
 7379       vBoundAtoms[*pat1]=set0;
 
 7380       for(vector<MolAtom*>::iterator pat2=this->GetAtomList().begin();pat2!=this->GetAtomList().end();++pat2)
 
 7384          for(list<StretchModeBondLength>::iterator pstretch=this->GetStretchModeBondLengthList().begin();
 
 7385             pstretch!=this->GetStretchModeBondLengthList().end();++pstretch)
 
 7387             set<MolAtom *>::iterator pos1=pstretch->mvTranslatedAtomList.find(*pat1),
 
 7388                                      pos2=pstretch->mvTranslatedAtomList.find(*pat2);
 
 7389             if(  ((pos1==pstretch->mvTranslatedAtomList.end())&&(pos2!=pstretch->mvTranslatedAtomList.end()))
 
 7390                ||((pos1!=pstretch->mvTranslatedAtomList.end())&&(pos2==pstretch->mvTranslatedAtomList.end())))
 
 7392                vBoundAtoms[*pat1].erase(*pat2);
 
 7401          for(list<StretchModeBondAngle>::iterator pstretch=this->GetStretchModeBondAngleList().begin();
 
 7402             pstretch!=this->GetStretchModeBondAngleList().end();++pstretch)
 
 7405             if((*pat1==pstretch->mpAtom1)||(*pat2==pstretch->mpAtom1)) 
continue;
 
 7407             set<MolAtom *>::iterator pos1=pstretch->mvRotatedAtomList.find(*pat1),
 
 7408                                      pos2=pstretch->mvRotatedAtomList.find(*pat2);
 
 7410             if(  ((pos1==pstretch->mvRotatedAtomList.end())&&(pos2!=pstretch->mvRotatedAtomList.end()))
 
 7411                ||((pos1!=pstretch->mvRotatedAtomList.end())&&(pos2==pstretch->mvRotatedAtomList.end())))
 
 7413                vBoundAtoms[*pat1].erase(*pat2);
 
 7422          for(list<StretchModeTorsion>::iterator pstretch=this->GetStretchModeTorsionList().begin();
 
 7423             pstretch!=this->GetStretchModeTorsionList().end();++pstretch)
 
 7425             set<MolAtom *>::iterator pos1=pstretch->mvRotatedAtomList.find(*pat1),
 
 7426                                      pos2=pstretch->mvRotatedAtomList.find(*pat2);
 
 7428             if(  (pos1!=pstretch->mvRotatedAtomList.end())&&(pos2==pstretch->mvRotatedAtomList.end())
 
 7429                &&(*pat2!=pstretch->mpAtom1) &&(*pat2!=pstretch->mpAtom2) )
 
 7431                vBoundAtoms[*pat1].erase(*pat2);
 
 7436             if(  (pos1==pstretch->mvRotatedAtomList.end())&&(pos2!=pstretch->mvRotatedAtomList.end())
 
 7437                &&(*pat1!=pstretch->mpAtom1) && (*pat1!=pstretch->mpAtom2) )
 
 7439                vBoundAtoms[*pat1].erase(*pat2);
 
 7449    set<set<MolAtom*> > vBoundGroups;
 
 7450    for(map<
MolAtom*,set<MolAtom*> >::iterator pos=vBoundAtoms.begin();pos!=vBoundAtoms.end();++pos)
 
 7453       cout<<
"Non-flexible group from "<<pos->first->GetName()<<
": ";
 
 7454       for(set<MolAtom*>::const_iterator atom=pos->second.begin();atom!=pos->second.end();++atom)
 
 7455          cout<<(*atom)->GetName()<<
",";
 
 7460          for(set<MolAtom *>::iterator at=(*pr)->begin();at!=(*pr)->end();++at)
 
 7461             pos->second.erase(*at);
 
 7462       if(pos->second.size()>1) vBoundGroups.insert(pos->second);
 
 7465    for(set<set<MolAtom*> >::iterator pos=vBoundGroups.begin();pos!=vBoundGroups.end();++pos)
 
 7467       cout<<
"Non-flexible group:";
 
 7468       for(set<MolAtom*>::const_iterator atom=pos->begin();atom!=pos->end();++atom)
 
 7469          cout<<(*atom)->GetName()<<
",";
 
 7475    for(set<set<MolAtom*> >::iterator pos=vBoundGroups.begin();pos!=vBoundGroups.end();++pos)
 
 7478       for(vector<MolBond*>::iterator pr=this->GetBondList().begin();pr!=this->GetBondList().end();++pr)
 
 7479          if(  (pos->find(&(*pr)->GetAtom1())!=pos->end())
 
 7480             ||(pos->find(&(*pr)->GetAtom2())!=pos->end())) vb.insert(*pr);
 
 7482       set<MolBondAngle*> va;
 
 7483       for(vector<MolBondAngle*>::iterator pr=this->GetBondAngleList().begin();pr!=this->GetBondAngleList().end();++pr)
 
 7484          if(  (pos->find(&(*pr)->GetAtom1())!=pos->end())
 
 7485             ||(pos->find(&(*pr)->GetAtom2())!=pos->end())
 
 7486             ||(pos->find(&(*pr)->GetAtom3())!=pos->end())) va.insert(*pr);
 
 7488       set<MolDihedralAngle*> vd;
 
 7489       for(vector<MolDihedralAngle*>::iterator pr=this->GetDihedralAngleList().begin();pr!=this->GetDihedralAngleList().end();++pr)
 
 7490          if(  (pos->find(&(*pr)->GetAtom1())!=pos->end())
 
 7491             ||(pos->find(&(*pr)->GetAtom2())!=pos->end())
 
 7492             ||(pos->find(&(*pr)->GetAtom3())!=pos->end())
 
 7493             ||(pos->find(&(*pr)->GetAtom4())!=pos->end())) vd.insert(*pr);
 
 7495       set<MolAtom*> tmp= *pos;
 
 7504    for(vector<MolAtom*>::iterator at=this->GetAtomList().begin();at!=this->GetAtomList().end();++at)
 
 7507       for(set<MolAtom *>::iterator at=(*pr)->begin();at!=(*pr)->end();++at)
 
 7509    cout<<
"Full MD atom group:"<<endl<<
" ";
 
 7511       cout<<(*pos)->GetName()<<
" ";
 
 7517       for(set<MolAtom*>::const_iterator at=pos->mvpAtom.begin();at!=pos->mvpAtom.end();++at)
 
 7519    cout<<
"Full MD atom group:"<<endl<<
" ";
 
 7521       cout<<(*pos)->GetName()<<
" ";
 
 7524    mClockMDAtomGroup.
Click();
 
 7534    VFN_DEBUG_ENTRY(
"Molecule::UpdateScattCompList()",5)
 
 7535    TAU_PROFILE(
"Molecule::UpdateScattCompList()",
"void ()",TAU_DEFAULT);
 
 7538    for(
long i=0;i<nb;++i)
 
 7548   #ifdef RIGID_BODY_STRICT_EXPERIMENTAL 
 7554         (*pos)->mQuat.Normalize();
 
 7556         REAL x0=0,y0=0,z0=0;
 
 7557         for(set<unsigned int>::iterator at=(*pos)->mvIdx.begin();at!=(*pos)->mvIdx.end();++at)
 
 7568         for(set<unsigned int>::iterator at=(*pos)->mvIdx.begin();at!=(*pos)->mvIdx.end();++at)
 
 7571           (*pos)->mQuat.RotateVector(x,y,z);
 
 7580    REAL x0=0,y0=0,z0=0;
 
 7583       for(
long i=0;i<nb;++i)
 
 7599    for(
long i=0;i<nb;++i)
 
 7608    for(
long i=0;i<nb;++i)
 
 7614    for(
long i=0;i<nb;++i)
 
 7621    for(
long i=0;i<nb;++i)
 
 7628    VFN_DEBUG_EXIT(
"Molecule::UpdateScattCompList()",5)
 
 7632    VFN_DEBUG_ENTRY(
"Molecule::FindAtom():"<<name,4)
 
 7633    vector<MolAtom*>::reverse_iterator rpos;
 
 7635       if(name==(*rpos)->GetName())
 
 7637          VFN_DEBUG_EXIT(
"Molecule::FindAtom():"<<name<<
"...NOT FOUND !",4)
 
 7640    VFN_DEBUG_EXIT(
"Molecule::FindAtom():"<<name<<
"...NOT FOUND !",4)
 
 7645    vector<MolAtom*>::const_reverse_iterator rpos;
 
 7648       if(name==(*rpos)->GetName()) 
return rpos;
 
 7653    VFN_DEBUG_ENTRY(
"Molecule::InitOptions",7)
 
 7654    static string Flexname;
 
 7655    static string Flexchoices[3];
 
 7657    static string FlipName;
 
 7658    static string FlipChoice[2];
 
 7660    static string autoOptimizeConformationName;
 
 7661    static string autoOptimizeConformationChoices[2];
 
 7663    static string optimizeOrientationName;
 
 7664    static string optimizeOrientationChoices[2];
 
 7666    static string moleculeCenterName;
 
 7667    static string moleculeCenterChoices[2];
 
 7669    static bool needInitNames=
true;
 
 7670    if(
true==needInitNames)
 
 7672       Flexname=
"Flexibility Model";
 
 7673       Flexchoices[0]=
"Automatic from Restraints, relaxed - RECOMMENDED";
 
 7674       Flexchoices[1]=
"Rigid Body";
 
 7675       Flexchoices[2]=
"Automatic from Restraints, strict";
 
 7678       FlipName=
"Enable Flipping";
 
 7679       FlipChoice[0]=
"Yes";
 
 7682       autoOptimizeConformationName=
"Auto Optimize Starting Conformation";
 
 7683       autoOptimizeConformationChoices[0]=
"Yes";
 
 7684       autoOptimizeConformationChoices[1]=
"No";
 
 7686       optimizeOrientationName=
"Optimize Orientation";
 
 7687       optimizeOrientationChoices[0]=
"Yes";
 
 7688       optimizeOrientationChoices[1]=
"No";
 
 7690       moleculeCenterName=
"Rotation Center";
 
 7691       moleculeCenterChoices[0]=
"Geometrical center (recommended)";
 
 7692       moleculeCenterChoices[1]=
"User-chosen Atom";
 
 7694       needInitNames=
false;
 
 7700    mFlipModel.Init(2, &FlipName, FlipChoice);
 
 7701    mFlipModel.SetChoice(0);
 
 7705                                   autoOptimizeConformationChoices);
 
 7714    VFN_DEBUG_EXIT(
"Molecule::InitOptions",7)
 
 7718 mpAtom0(&at0),mpAtom1(&at1),mpAtom2(&at2),mNbTest(0),mNbAccept(0)
 
 7724    TAU_PROFILE(
"Molecule::FlipAtomGroup(FlipGroup&)",
"void (...)",TAU_DEFAULT);
 
 7738       const REAL norm01=sqrt(v01x*v01x+v01y*v01y+v01z*v01z+1e-7);
 
 7739       v01x /= norm01;v01y /= norm01;v01z /= norm01;
 
 7744       const REAL norm02=sqrt(v02x*v02x+v02y*v02y+v02z*v02z+1e-7);
 
 7745       v02x /= norm02;v02y /= norm02;v02z /= norm02;
 
 7750       const REAL norm12=sqrt(v12x*v12x+v12y*v12y+v12z*v12z+1e-7);
 
 7751       v12x /= norm12;v12y /= norm12;v12z /= norm12;
 
 7756       const REAL norm0m=sqrt(v0mx*v0mx+v0my*v0my+v0mz*v0mz+1e-7);
 
 7757       v0mx /= norm0m;v0my /= norm0m;v0mz /= norm0m;
 
 7759       if(fabs(v01x*v02x+v01y*v02y+v01z*v02z)
 
 7760          >0.05*sqrt(abs( (v01x*v01x+v01y*v01y+v01z*v01z)
 
 7761                         *(v02x*v02x+v02y*v02y+v02z*v02z))))
 
 7763          REAL v012x=v01y*v02z-v01z*v02y;
 
 7764          REAL v012y=v01z*v02x-v01x*v02z;
 
 7765          REAL v012z=v01x*v02y-v01y*v02x;
 
 7766          const REAL norm012=sqrt(v012x*v012x+v012y*v012y+v012z*v012z+1e-7);
 
 7767          v012x /= norm012;v012y /= norm012;v012z /= norm012;
 
 7770          for(list<pair<
const MolAtom *,set<MolAtom*> > >::const_iterator
 
 7774             REAL v03x=chain->first->X()-group.
mpAtom0->X();
 
 7775             REAL v03y=chain->first->Y()-group.
mpAtom0->Y();
 
 7776             REAL v03z=chain->first->Z()-group.
mpAtom0->Z();
 
 7777             const REAL norm03=sqrt( v03x*v03x + v03y*v03y + v03z*v03z +1e-7);
 
 7778             v03x /= norm03;v03y /= norm03;v03z /= norm03;
 
 7780             const REAL a1=v012x*v03x+v012y*v03y+v012z*v03z;
 
 7781             const REAL a2= v0mx*v03x+ v0my*v03y+ v0mz*v03z;
 
 7782             const REAL a3= v12x*v03x+ v12y*v03y+ v12z*v03z;
 
 7783             REAL angle = -a1/sqrt(1-a3*a3+1e-7);
 
 7792                else angle = asin(angle);
 
 7794             if(a2<0) angle=M_PI-angle;
 
 7796                                   chain->second,2*angle,keepCenter);
 
 7808    #ifdef RIGID_BODY_STRICT_EXPERIMENTAL 
 7813       (*pos)->mQuat.Normalize();
 
 7815       REAL x0=0,y0=0,z0=0;
 
 7816       for(set<MolAtom *>::iterator at=(*pos)->begin();at!=(*pos)->end();++at)
 
 7827       for(set<MolAtom *>::iterator at=(*pos)->begin();at!=(*pos)->end();++at)
 
 7829         REAL x=(*at)->GetX()-x0, y=(*at)->GetY()-y0, z=(*at)->GetZ()-z0;
 
 7830         (*pos)->mQuat.RotateVector(x,y,z);
 
 7831         (*at)->SetX(x+x0+(*pos)->mX);
 
 7832         (*at)->SetY(y+y0+(*pos)->mY);
 
 7833         (*at)->SetZ(z+z0+(*pos)->mZ);
 
 7840       (*pos)->mQuat.Q0()=1;
 
 7841       (*pos)->mQuat.Q1()=0;
 
 7842       (*pos)->mQuat.Q2()=0;
 
 7843       (*pos)->mQuat.Q3()=0;
 
 7848 #ifdef __WX__CRYST__ 
 7851    VFN_DEBUG_ENTRY(
"Molecule::WXCreate()",5)
 
 7853    VFN_DEBUG_EXIT("
Molecule::WXCreate()",5)
 
 7854    return mpWXCrystObj;
 
MolAtom(const REAL x, const REAL y, const REAL z, const ScatteringPower *pPow, const string &name, Molecule &parent)
Constructor for a MolAtom. 
CrystVector_REAL GetLatticePar() const 
Lattice parameters (a,b,c,alpha,beta,gamma) as a 6-element vector in Angstroems and radians...
RotorGroup(const MolAtom &at1, const MolAtom &at2)
Constructor, with the two atoms around which the rotation shall be made. 
void SetDerivStep(const REAL)
Fixed step to use to compute numerical derivative. 
list< MDAtomGroup > mvMDAtomGroup
Groups of atoms that should be moved according to molecular dynamics principles. 
REAL mOccupancy
Occupancy. 
void RestraintExport(ostream &os) const 
Print the restraints (bond length, angles...) as whole labels and number in column text format which ...
map< MolAtom *, set< MolAtom * > > mConnectivityTable
Connectivity table: for each atom, keep the list of atoms bonded to it. 
MolAtom * mpAtom1
The second atom (first atom moved) 
StretchModeTwist(MolAtom &at1, MolAtom &at2)
Constructor If pDihedralAngle!=0, the dihedral angle length restraint is respected. 
XYZ mDerivAtom1
Partial derivatives of the angle with respect to the coordinates of the atoms. 
void AddPar(const RefinablePar &newRefPar)
Add a refinable parameter. 
vector< MolDihedralAngle * >::iterator RemoveDihedralAngle(const MolDihedralAngle &, const bool del=true)
Remove a dihedral angle. 
virtual void UpdateDisplay() const 
If there is an interface, this should be automatically be called each time there is a 'new...
vector< MolBond * > mvpBond
The list of bonds. 
REAL mXmin
Display limits in reduced coordinates. 
bool mShowLabel
Show labels ? 
void BuildRotorGroup()
Build the groups of atoms that will be rotated during global optimization. 
REAL mLLK
Stored log(likelihood) 
REAL mX
The translation of all the atoms as a group The values will be resetted whenever entering or leaving ...
Bond angle restraint between 3 atoms. 
void TuneGlobalOptimRotationAmplitude()
Tune the rotation amplitude for free torsions and for the overall Molecule Rotation. 
virtual void Print(ostream &os, bool full=true) const 
Print one-line list of atoms moved. 
vector< MolDihedralAngle * >::const_iterator FindDihedralAngle(const MolAtom &at1, const MolAtom &at2, const MolAtom &at3, const MolAtom &at4) const 
Searches whether a dihedral between four atoms already exists, searching for either (at1...
Molecule * mpMol
Parent Molecule. 
virtual void Stretch(const REAL change, const bool keepCenter=true)
Move the atoms according to this mode. 
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. 
virtual REAL GetLogLikelihood() const 
Get -ln(likelihood) for this restraint. 
float string2floatC(const string &s)
Function to convert a substring to a floating point value, imposing a C locale (using '...
bool mIsInRing
Is the atom in a ring ? 
Atoms moved when changing a bond angle. 
REAL GetBondLength(const MolAtom &at1, const MolAtom &at2)
Get The Bond Length between two atoms. 
MolBondAngle(MolAtom &atom1, MolAtom &atom2, MolAtom &atom3, const REAL angle, const REAL sigma, const REAL delta, Molecule &parent)
Constructor. 
void OrthonormalToFractionalCoords(REAL &x, REAL &y, REAL &z) const 
Get fractional cartesian coordinates for a set of (x,y,z) orthonormal coordinates. 
virtual void CalcDeriv(const bool derivllk=true) const 
Calculate the derivative of the Molecule's Log(likelihood) and atomic positions versus a change of th...
virtual const string & GetClassName() const 
Name for this class ("RefinableObj", "Crystal",...). 
set< MolAtom * > mvRotatedAtomList
The set of atoms that are to be rotated around the direction going through at1 and perpendicular to t...
long GetNbComponent() const 
Number of components. 
We need to record exactly when refinable objects have been modified for the last time (to avoid re-co...
std::map< const MolAtom *, XYZ > mDerivXYZ
Derivative of the atomic positions versus a change of the bond length. 
virtual ~MolAtom()
Destructor. 
bool IsBeingRefined() const 
Is the object being refined ? (Can be refined by one algorithm at a time only.) 
void BuildFlipGroup()
Build the groups of atoms that can be flipped. 
virtual void GlobalOptRandomMove(const REAL mutationAmplitude, const RefParType *type=gpRefParTypeObjCryst)
Make a random move of the current configuration. 
void FractionalToOrthonormalCoords(REAL &x, REAL &y, REAL &z) const 
Get orthonormal cartesian coordinates for a set of (x,y,z) fractional coordinates. 
string GetName() const 
Name of the bond, e.g. "C3-O4". 
CrystVector_REAL mXYZ
coordinates of the scatterer (or of its center..) 
vector< MolBond * >::const_iterator FindBond(const MolAtom &, const MolAtom &) const 
Searches whether a bond between two atoms already exists. 
CrystVector_REAL mLSQObs
Current LSQ Calc - one value for each restraint (bond distance, angle or dihedral angle ideal values)...
void AddRefinableObj(RefinableObj &)
Add a refined object. All sub-objects are also added. 
MolAtom * mpAtom2
The third atom. 
std::vector< RigidGroup * >::iterator RemoveRigidGroup(const RigidGroup &group, const bool updateDisplay=true, const bool del=true)
Remove a rigid group of atoms. 
REAL mLLK
Stored log(likelihood) 
virtual void Stretch(const REAL change, const bool keepCenter=true)
Move the atoms according to this mode. 
void BuildStretchModeTorsion()
Build the groups of atoms moved when changing a dihedral angle, while respecting the Molecule restrai...
REAL mDerivLLKCoeff
The factor used to change the derivative of the length/angle, to the derivative of the log(likelihood...
virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false)
This should be called by any optimization class at the begining of an optimization. 
void SaveParamSet(const unsigned long id) const 
Save the current set of refined values over a previously-created set of saved values. 
const MolAtom * mpCenterAtom
Atom chosen as center of rotation, if mRotationCenter is set to use an atom rather than the geometric...
StretchModeBondLength(MolAtom &at0, MolAtom &at1, const MolBond *pBond)
Constructor If pBond!=0, the bond length restraint is respected. 
list< StretchModeBondLength > mvStretchModeBondLength
List of StretchModeBondLength. 
virtual void GlobalOptRandomMove(const REAL mutationAmplitude, const RefParType *type)
Make a random move of the current configuration. 
vector< MolAtom * > mvpAtom
The vector of the 3 atoms involved in the bond angle. 
Molecule * mpMol
The Molecule corresponding to this stretch mode. 
virtual void TagNewBestConfig() const 
During a global optimization, tells the object that the current config is the latest "best" config...
MolAtom * mpAtom1
The first atom. 
const MolAtom * mpAtom1
The first atom defining the rotation axis. 
virtual void Print(ostream &os, bool full=true) const 
Print one-line list of atoms moved. 
void AddBondAngle(MolAtom &atom1, MolAtom &atom2, MolAtom &atom3, const REAL angle, const REAL sigma, const REAL delta, const bool updateDisplay=true)
Add a bond angle restraint. 
virtual REAL GetLogLikelihood() const 
Get -ln(likelihood) for this restraint. 
void Click()
Record an event for this clock (generally, the 'time' an object has been modified, or some computation has been made) 
FlipGroup(const MolAtom &at0, const MolAtom &at1, const MolAtom &at2)
Constructor, with the central atom. 
vector< Restraint * >::iterator RemoveRestraint(Restraint *pRestraint)
Remove a restraint from the list of known restraints. 
void AddChild(const RefinableObjClock &)
Add a 'child' clock. 
bool mDeleteSubObjInDestructor
Base Rotation amplitude (in radians) for the Molecule, so that the average atomic displacement is equ...
virtual const CrystVector_REAL & GetLSQDeriv(const unsigned int n, RefinablePar &par)
Get the first derivative values for the LSQ function, for a given parameter. 
Quaternion mQuat
The unit quaternion defining the orientation - this is used during optimizations to rotate all atoms ...
void ResetRigidGroupsPar() const 
Set the orientation & translation parameters of all rigid groups to 0, after correcting the atomic po...
virtual const CrystVector_REAL & GetLSQObs(const unsigned int) const 
Get the observed values for the LSQ function. 
REAL DihedralAngleRandomChange(const StretchModeTorsion &mode, const REAL amplitude, const bool respectRestraint=true)
Change a dihedral angle, while respecting the Restraint (if any). 
vector< MolAtom * > mvpAtom
The list of atoms. 
void XMLInput(istream &is, const XMLCrystTag &tag)
XMLInput From stream. 
std::vector< RigidGroup * > mvRigidGroup
Rigid groups of atoms. 
RefinablePar & GetPar(const long i)
Access all parameters in the order they were inputted. 
StretchModeBondAngle(MolAtom &at0, MolAtom &at1, MolAtom &at2, const MolBondAngle *pBondAngle)
Constructor If pBondAngle!=0, the bond angle length restraint is respected. 
A quaternion class, used to represent the orientation of the molecule. 
const MolAtom * mpAtom0
The atom which is an asymmetric center. 
Bond between two atoms, also a restraint on the associated bond length. 
REAL mLLK
Stored log(likelihood) 
virtual ostream & POVRayDescription(ostream &os, const CrystalPOVRayOptions &options) const 
ScatteringComponentList mScattCompList
The list of scattering components. 
MDAtomGroup()
Default constructor. 
void Print(ostream &os, bool full=true) const 
Print one-line list of atoms moved. 
std::list< StretchMode * > mvpStretchModeNotFree
Groups of StretchMode breaking restraints (beyond the one they are associated to) ...
wxCryst class for Molecule objects 
virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false)
This should be called by any optimization class at the begining of an optimization. 
void RotateAtomGroup(const MolAtom &at1, const MolAtom &at2, const set< MolAtom * > &atoms, const REAL angle, const bool keepCenter=true)
Rotate a group of atoms around an axis defined by two atoms. 
MolDihedralAngle(MolAtom &atom1, MolAtom &atom2, MolAtom &atom3, MolAtom &atom4, const REAL angle, const REAL sigma, const REAL delta, Molecule &parent)
Constructor. 
list< MolRing > mvRing
The list of rings. 
Structure holding 3 coordinates, or deriviatives with respect to each of these coordinates. 
XYZ mDerivAtom1
Derivatives of the bond length with respect to the coordinates of the atoms. 
Atoms moved when rotated around a bond at0-at1-at2-at3. 
MolAtom * mpAtom2
The second atom. 
list< StretchModeTwist > mvStretchModeTwist
List of StretchModeTwist. 
void AddDihedralAngle(MolAtom &atom1, MolAtom &atom2, MolAtom &atom3, MolAtom &atom4, const REAL angle, const REAL sigma, const REAL delta, const bool updateDisplay=true)
Add a dihedral angle restraint. 
Rigid groups of atoms inside a molecule. 
virtual int GetNbComponent() const 
Number of components in the scatterer (eg number of point scatterers) 
RefinableObjClock & GetBondListClock()
get the clock associated to the list of bonds 
void AddAtom(const REAL x, const REAL y, const REAL z, const ScatteringPower *pPow, const string &name, const bool updateDisplay=true)
Add an atom. 
unsigned int GetNbOption() const 
Number of Options for this object. 
RefObjOpt & GetOption(const unsigned int i)
Access to the options. 
std::vector< MolZAtom > mAsZMatrix
The Molecule, as a lightweight ZMatrix, for export purposes. 
void CalcGradient(std::map< MolAtom *, XYZ > &m) const 
Calc log(likelihood) gradient - versus all atomic coordinates. 
list< StretchModeTorsion > mvStretchModeTorsion
List of StretchModeBondLength. 
REAL mQ0
The components of the quaternion z=(q0,v) with v=(q1,q2,q3) 
list< RotorGroup > mvRotorGroupTorsion
List of RotorGroups corresponding to free torsion bonds. 
CrystVector_REAL mLSQCalc
Current LSQ Calc - one value for each restraint (bond distance, angle or dihedral angle) ...
MolBond(MolAtom &atom1, MolAtom &atom2, const REAL length, const REAL sigma, const REAL delta, Molecule &parent, const REAL bondOrder=1.)
Constructor. 
void BuildStretchModeTwist()
Build the groups of atoms used to twist internally the Molecule, e.g. 
virtual void Print(ostream &os, bool full=true) const 
Print one-line list of atoms moved. 
void MolecularDynamicsEvolve(std::map< MolAtom *, XYZ > &v0, const unsigned nbStep, const REAL dt, const std::vector< MolBond * > &vb, const std::vector< MolBondAngle * > &va, const std::vector< MolDihedralAngle * > &vd, std::map< RigidGroup *, std::pair< XYZ, XYZ > > &vr, REAL nrj0=0)
Change the conformation of the molecule using molecular dynamics principles. 
REAL GetZ() const 
Z coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
std::map< const MolDihedralAngle *, REAL > mvpBrokenDihedralAngle
List of dihedral angle restraints modified by this mode The key is the restraint, the value is the de...
void BuildStretchModeBondAngle()
Build the groups of atoms moved when changing a bond angle, while respecting the Molecule restraints...
void BuildRingRecursive(MolAtom *currentAtom, MolAtom *previousAtom, const map< MolAtom *, set< MolAtom * > > &connect, list< MolAtom * > &atomlist, map< set< MolAtom * >, list< MolAtom * > > &ringlist)
Find rings, starting from a one atom, and given a connectivity table. 
ObjRegistry< Crystal > gCrystalRegistry("List of all Crystals")
Global registry for all Crystal objects. 
REAL GetOccupancy() const 
Get the occupancy of the scatterer (0. 
const MolAtom * GetCenterAtom() const 
Get the atom defining the origin of the Molecule Equal to 0 if no atom as been set. 
Abstract base class for all objects in wxCryst. 
virtual const CrystVector_REAL & GetLSQWeight(const unsigned int) const 
Get the weight values for the LSQ function. 
Base object for Monte-Carlo Global Optimization methods. 
Quaternion operator*(const Quaternion &q) const 
Quaternion multiplication. 
const ScatteringPower * mpScattPow
ScatteringPower. 
virtual void RandomStretch(const REAL amplitude, const bool keepCenter=true)
Move the atoms according to this mode, randomly. 
string mName
Name for this RefinableObject. Should be unique, at least in the same scope.+. 
void BuildStretchModeGroups()
Separate StretchMode that break more than their assigned restraint from others. 
const std::vector< MolZAtom > & AsZMatrix(const bool keeporder) const 
Molecule as Z-matrix. 
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 
virtual const CrystVector_REAL & GetLSQDeriv(const unsigned int, RefinablePar &)
Get the first derivative values for the LSQ function, for a given parameter. 
void BuildMDAtomGroups()
Find groups of atoms that cannot be moved relatively to each other using the free or non-free stretch...
virtual Molecule * CreateCopy() const 
static Quaternion RotationQuaternion(const REAL ang, const REAL v1, const REAL v2, const REAL v3)
Create a rotation quaternion around a given vector for a given angle. 
Quaternion GetConjugate() const 
Get the conjugate of this quaternion (== the inverse if unit quaternion) 
void SetIsInRing(const bool r) const 
Flag this atom as being in a ring (or not). 
virtual REAL GetLogLikelihood() const 
Get -ln(likelihood) for this restraint. 
virtual void RandomStretch(const REAL amplitude, const bool keepCenter=true)
Move the atoms according to this mode, randomly. 
virtual ~MolBondAngle()
Destructor. 
REAL mX
Cartesian oordinates in the Molecule reference frame. 
REAL mDerivLLKCoeff
The factor used to change the derivative of the length/angle, to the derivative of the log(likelihood...
REAL mLogLikelihood
The current log(likelihood) 
REAL GetBondAngle(const MolAtom &at1, const MolAtom &at2, const MolAtom &at3)
Get The Bond Angle of 3 atoms. 
virtual void Stretch(const REAL change, const bool keepCenter=true)
Move the atoms according to this mode. 
Defines a group of atoms which can be rotated around an axis defined by two other atoms...
Molecule : class for complex scatterer descriptions using cartesian coordinates with bond length/angl...
vector< MolAtom * >::reverse_iterator FindAtom(const string &name)
Search a MolAtom from its name. 
Class to store POV-Ray output options. 
Molecule(Crystal &cryst, const string &name="")
Constructor. 
Molecule * mpMol
Parent Molecule. 
const MolDihedralAngle * mpDihedralAngle
The (optional) bond angle restraint which this stretch mode should respect. 
void BuildStretchModeBondLength()
Build the groups of atoms moved when stretching a bond length, while respecting the Molecule restrain...
void RestoreParamSet(const unsigned long id)
Restore a saved set of values. 
vector< MolBondAngle * >::iterator RemoveBondAngle(const MolBondAngle &, const bool del=true)
Remove a BondAngle. 
set< MolAtom * > mvTranslatedAtomList
The set of atoms that are to be translated, including at1. 
void FlipAtomGroup(const FlipGroup &, const bool keepCenter=true)
Flip a group of atom. See Molecule::FlipGroup. 
void SetAlgorithmParallTempering(const AnnealingSchedule scheduleTemp, const REAL tMax, const REAL tMin, const AnnealingSchedule scheduleMutation=ANNEALING_CONSTANT, const REAL mutMax=16., const REAL mutMin=.125)
Set the refinement method to Parallel Tempering. 
RefinableObjClock mClockScattCompList
void AssignClock(RefinableObjClock &clock)
bool mShowHydrogens
Show hydrogens ? 
vector< Restraint * > mvpRestraint
Vector of pointers to the restraints for this object. 
virtual string GetComponentName(const int i) const 
Name for the i-th component of this scatterer. 
void XMLOutput(ostream &os, int indent=0) const 
XMLOutput to stream in well-formed XML. 
REAL GetX() const 
X coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
Vector library (Blitz++ mimic) for ObjCryst++. 
void UpdateScattCompList() const 
Update the Molecule::mScattCompList from the cartesian coordinates of all atoms, and the orientation ...
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input From stream. 
virtual unsigned int GetNbLSQFunction() const 
Number of LSQ functions. 
REAL mDerivLLKCoeff
The factor used to change the derivative of the length/angle, to the derivative of the log(likelihood...
Molecule * mpMol
Parent Molecule. 
MolAtom * mpAtom2
The second atom. 
void Normalize() const 
Re-normalize the quaternion to unity. 
REAL mOccupancy
Occupancy : 0 <= occ <= 1 For a multi-atom scatterer (polyhedron,..), this is the overall occupancy o...
XYZ mDerivAtom1
Partial derivatives of the angle with respect to the coordinates of the atoms. 
REAL LorentzianBiasedRandomMove(const REAL x0, const REAL sigma, const REAL delta, const REAL amplitude)
Random move respecting a gaussian probability distribution with a flat top. 
void CalcGradient(std::map< MolAtom *, XYZ > &m) const 
Calc log(likelihood) gradient - versus all atomic coordinates. 
REAL GetDynPopCorr(const Scatterer *pscatt, unsigned int component) const 
Access the Dynamical Occupancy Correction for a given component (atom) in a given Scatterer...
void RotateVector(REAL &v1, REAL &v2, REAL &v3) const 
Rotate vector v=(v1,v2,v3). The rotated components are directly written. 
virtual void RandomizeConfiguration()
Randomize Configuration (before a global optimization). 
MolAtom : atom inside a Molecule. 
void SetCenterAtom(const MolAtom &at)
Get the atom defining the origin of the Molecule Equal to 0 if no atom as been set. 
REAL mLLKDeriv
Derivative of the Molecule's Log(likelihood) versus a change of the bond length. 
REAL GetDihedralAngle(const MolAtom &at1, const MolAtom &at2, const MolAtom &at3, const MolAtom &at4)
Get The dihedral angle defined by 4 atoms. 
REAL GetDeriv(const std::map< const MolAtom *, XYZ > &m, const bool llk=false) const 
Get the derivative of the bond length, given the derivatives of the atom positions This requires that...
REAL mBaseAmplitude
The recommended change amplitude, for a base global optimization displacement, to obtain an average 0...
const map< MolAtom *, set< MolAtom * > > & GetConnectivityTable()
Get the connectivity table. 
RefObjOpt mFlexModel
OPtion for the different types of flexibility possible for this molecule: rigid body, free atoms + restraints, torsion angles... 
void InitOptions()
Build options for this object. 
vector< MolAtom * >::iterator RemoveAtom(MolAtom &, const bool del=true)
Remove an atom. 
virtual ~MolDihedralAngle()
Destructor. 
REAL mLogLikelihoodScale
Scale (multiplier) for the log(likelihood) 
virtual void InitRefParList()
Quaternion mQuat
The unit quaternion defining the orientation. 
void AddRigidGroup(const RigidGroup &, const bool updateDisplay=true)
Add a rigid group of atoms. 
MolAtom * mpAtom0
The first atom. 
const MolBond * mpBond
The (optional) bond length which this stretch mode should respect. 
virtual void XMLOutput(ostream &os, int indent=0) const 
Output to stream in well-formed XML. 
REAL GetDeriv(const std::map< const MolAtom *, XYZ > &m, const bool llk=false) const 
Get the derivative of the angle, given the derivatives of the atom positions This requires that GetLo...
RefinableObjClock & GetRigidGroupClock()
Get the clock associated to the list of rigid groups (clicked also whenever a rigid group is modified...
list< RotorGroup > mvRotorGroupInternal
List of RotorGroups for internal rotations. 
const MolBondAngle * mpBondAngle
The (optional) bond angle restraint which this stretch mode should respect. 
Molecule * mpMol
Parent Molecule. 
string mName
Name for this atom. 
void SetDeleteSubObjInDestructor(const bool b)
Set whether to delete the MolAtoms, MolBonds, MolBondAngles and MolDihedralAngles in the destructor...
void CalcGradient(std::map< MolAtom *, XYZ > &m) const 
Calc log(likelihood) gradient - versus all atomic coordinates. 
std::map< const MolBond *, REAL > mvpBrokenBond
List of bond restraints affected by this mode The key is the restraint, the value is the derivative o...
list< RotorGroup > mvRotorGroupTorsionSingleChain
List of RotorGroups corresponding to free torsion bonds, but with only one chain of atoms listed...
Dihedral angle restraint between 4 atoms. 
void ClearParamSet(const unsigned long id) const 
Erase the param set with the given id, releasing memory. 
virtual void RandomizeConfiguration()
Randomize Configuration (before a global optimization). 
vector< MolBondAngle * > mvpBondAngle
The list of bond angles. 
virtual void SetName(const string &name)
Name of the object. 
std::list< StretchMode * > mvpStretchModeFree
Groups of StretchMode not breaking any restraint (unless the one they are associated to) ...
list< FlipGroup > mvFlipGroup
The list of FlipGroups. 
virtual void CalcDeriv(const bool derivllk=true) const 
Calculate the derivative of the Molecule's Log(likelihood) and atomic positions versus a change of th...
const Crystal & GetCrystal() const 
In which crystal is this Scatterer included ? 
virtual REAL GetLogLikelihood() const 
Get -log(likelihood) of the current configuration for the object. 
bool IsDescendantFromOrSameAs(const RefParType *type) const 
Returns true if the parameter is a descendant of 'type'. 
void Reset()
Reset a Clock to 0, to force an update. 
virtual void RandomStretch(const REAL amplitude, const bool keepCenter=true)
Move the atoms according to this mode, randomly. 
Exception class for ObjCryst++ library. 
vector< MolBondAngle * >::const_iterator FindBondAngle(const MolAtom &at1, const MolAtom &at0, const MolAtom &at2) const 
Searches whether a bond between three atoms already exists, searching for either (at1,at2,at3) and (at3,at2,at1), as these are equivalent. 
virtual void EndOptimization()
This should be called by any optimization class at the end of an optimization. 
virtual const CrystVector_REAL & GetLSQCalc(const unsigned int) const 
Get the current calculated value for the LSQ function. 
Cubic spline interpolation. 
RefObjOpt mAutoOptimizeConformation
Option to automatically optimize the starting conformation, if the total restraint cost is too high...
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++. 
REAL mMDMoveFreq
Frequency of using molecular dynamics move during GlobalOptRandomMove() 
virtual void EndOptimization()
This should be called by any optimization class at the end of an optimization. 
StretchModeTorsion(MolAtom &at1, MolAtom &at2, const MolDihedralAngle *pDihedralAngle)
Constructor If pDihedralAngle!=0, the dihedral angle length restraint is respected. 
std::map< const MolBondAngle *, REAL > mvpBrokenBondAngle
List of bond angle restraints modified by this mode The key is the restraint, the value is the deriva...
void BuildConnectivityTable() const 
Build the Connectivity table. 
std::string GetFormula() const 
Formula with atoms in alphabetic order. 
virtual REAL GetLogLikelihood() const 
Get -log(likelihood) of the current configuration for the object. 
Generic class for parameters of refinable objects. 
virtual const ScatteringComponentList & GetScatteringComponentList() const 
Get the list of all scattering components for this scatterer. 
REAL BondLengthRandomChange(const StretchModeBondLength &mode, const REAL amplitude, const bool respectRestraint=true)
Stretch a bond, while respecting the Restraint (if any). 
virtual const string & GetName() const 
Name of the object. 
virtual void Stretch(const REAL change, const bool keepCenter=true)
Move the atoms according to this mode. 
const std::vector< RigidGroup * > & GetRigidGroupList() const 
List of rigid group of atoms. 
virtual void RandomStretch(const REAL amplitude, const bool keepCenter=true)
Move the atoms according to this mode, randomly. 
set< MolAtom * > mvRotatedAtomList
The set of atoms that are to be rotated around at1-at2. 
void ExpandAtomGroupRecursive(MolAtom *atom, const map< MolAtom *, set< MolAtom * > > &connect, set< MolAtom * > &atomlist, const MolAtom *finalAtom)
Build recursively a list of atoms, starting from a one atom, and given a connectivity table...
void BuildRingList()
Build the list of rings in the molecule. 
virtual void Print(ostream &os, bool full=true) const 
Print one-line list of atoms moved. 
void RestraintStatus(ostream &os) const 
Print the status of all restraints (bond length, angles...) 
set< MolAtom * > mvRotatedAtomList
The set of atoms that are to be rotated around at1-at2. 
void AddRestraint(Restraint *pNewRestraint)
Add a new restraint. 
virtual ~MolBond()
Destructor. 
MolAtom * mpAtom1
The second atom. 
void SetX(const REAL) const 
Set the X,Y,Z coordinate - this is const because sometimes their coordinate must be changed even thou...
virtual void UpdateDisplay() const 
If there is an interface, this should be automatically be called each time there is a 'new...
RefinableObjClock & GetAtomPositionClock()
Get the clock associated to the atomic positions. 
Groups of atoms that can be moved using molecular dynamics principles, taking a list of restraints as...
RefObjOpt mMoleculeCenter
Option to choose the center of rotation of the Molecule for the global orientation either as the geom...
Crystal class: Unit cell, spacegroup, scatterers. 
vector< MolBond * >::iterator RemoveBond(const MolBond &, const bool del=true)
Remove a bond. 
vector< RefinablePar * >::iterator RemovePar(RefinablePar *refPar)
Remove a refinable parameter. 
vector< MolDihedralAngle * > mvpDihedralAngle
The list of dihedral angles. 
const MolAtom * mpAtom2
The second atom defining the rotation axis. 
REAL GetY() const 
Y coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
virtual void CalcDeriv(const bool derivllk=true) const 
Calculate the derivative of the Molecule's Log(likelihood) and atomic positions versus a change of th...
class to input or output a well-formatted xml beginning or ending tag. 
ObjRegistry< RefObjOpt > mOptionRegistry
List of options for this object. 
Group of atoms for random moves changing a bond length. 
CrystVector_REAL mLSQWeight
Current LSQ Calc - one value for each restraint(bond distance, angle or dihedral angle sigmas) ...
MolAtom * mpAtom1
The first atom. 
RefinableObjClock mClockScatterer
Last time anything (number of atoms, positions, scattering power) was changed. 
list< pair< const MolAtom *, set< MolAtom * > > > mvRotatedChainList
The set of atoms that are to be rotated during the flip. 
ScatteringPower & GetScatteringPower(const string &name)
Find a ScatteringPower from its name. Names must be unique in a given Crystal. 
class of refinable parameter types. 
list of scattering positions in a crystal, associated with the corresponding occupancy and a pointer ...
REAL GetDeriv(const std::map< const MolAtom *, XYZ > &m, const bool llk=false) const 
Get the derivative of the Angle, given the derivatives of the atom positions This requires that GetLo...
virtual void Print() const 
Print some info about the scatterer (ideally this should be one line...). 
unsigned long CreateParamSet(const string name="") const 
Save the current set of refined values in a new set. 
REAL mMDMoveEnergy
Relative energy of molecule during molecular dynamics move Default: 40, 10 (slow conformation change)...
void SetName(const string &)
Set the name of the parameter. It should be unique in the RefinableObj. 
virtual void CalcDeriv(const bool derivllk=true) const 
Calculate the derivative of the Molecule's Log(likelihood) and atomic positions versus a change of th...
Crystal * mpCryst
The crystal in which the Scatterer is This is needed so that we can know which scattering powers are ...
MolAtom * mpAtom0
The first atom (fixed). 
std::set< MolAtom * > mvMDFullAtomGroup
Full list of atoms that can be moved using molecular dynamics This excludes any atom part of a rigid ...
void RigidifyWithDihedralAngles()
Add dihedral angles so as to rigidify the Molecule. 
Quaternion()
Default constructor, yields q=(1,0,0,0) 
REAL BondAngleRandomChange(const StretchModeBondAngle &mode, const REAL amplitude, const bool respectRestraint=true)
change a bond angle, while respecting the Restraint (if any). 
void AddOption(RefObjOpt *opt)
void SetGlobalOptimStep(const REAL)
Maximum step to use during Global Optimization algorithms. 
list< StretchModeBondAngle > mvStretchModeBondAngle
List of StretchModeBondLength. 
void TranslateAtomGroup(const set< MolAtom * > &atoms, const REAL dx, const REAL dy, const REAL dz, const bool keepCenter=true)
Translate a group of atoms in a given direction. 
virtual void SetName(const string &name)
Name of the object. 
void AddBond(MolAtom &atom1, MolAtom &atom2, const REAL length, const REAL sigma, const REAL delta, const REAL bondOrder=1., const bool updateDisplay=true)
Add a bond. 
bool IsDummy() const 
Returns true if this is a dummy atom, i.e. 
Atoms moved between two other atoms, using a "twist" of their positions - only small twists of their ...
When 3(A1..1n) or more atoms are connected to a same atom A, it defines a 'flip' group, where it is possible to rotate bonds to their symmetric with respect to one plane defined by atoms Ai-A-Aj. 
void XMLOutput(ostream &os, const string &name, int indent=0) const 
XMLOutput to stream in well-formed XML. 
RefObjOpt mOptimizeOrientation
Option to optimize the Molecule's orientation. 
vector< MolAtom * > mvpAtom
The vector of the 4 atoms involved in the bond angle. 
void OptimizeConformationSteepestDescent(const REAL maxStep=0.1, const unsigned nbStep=1)
Optimize the conformation from internal restraints (bond lengths, angles and dihedral angles)...
void OptimizeConformation(const long nbTrial=10000, const REAL stopCost=0.)
Minimize configuration from internal restraints (bond lengths, angles and dihedral angles)...
Abstract Base Class to describe the scattering power of any Scatterer component in a crystal...