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...