28 #include <boost/format.hpp>
30 #include "ObjCryst/Quirks/VFNStreamFormat.h"
31 #include "ObjCryst/ObjCryst/Molecule.h"
32 #include "ObjCryst/ObjCryst/ZScatterer.h"
33 #include "ObjCryst/RefinableObj/GlobalOptimObj.h"
37 #include <OpenGL/glu.h>
44 #include "ObjCryst/wxCryst/wxMolecule.h"
50 #define RIGID_BODY_STRICT_EXPERIMENTAL
53 #undef RESTRAINT_X2_X4_X6
60 XYZ::XYZ(REAL x0,REAL y0,REAL z0):x(x0),y(y0),z(z0){};
91 return sqrt( (at1.GetX()-at2.GetX())
92 *(at1.GetX()-at2.GetX())
93 +(at1.GetY()-at2.GetY())
94 *(at1.GetY()-at2.GetY())
95 +(at1.GetZ()-at2.GetZ())
96 *(at1.GetZ()-at2.GetZ()) );
101 const REAL x21=at1.GetX()-at2.GetX();
102 const REAL y21=at1.GetY()-at2.GetY();
103 const REAL z21=at1.GetZ()-at2.GetZ();
104 const REAL x23=at3.GetX()-at2.GetX();
105 const REAL y23=at3.GetY()-at2.GetY();
106 const REAL z23=at3.GetZ()-at2.GetZ();
107 const REAL norm21_norm23= sqrt( (x21*x21+y21*y21+z21*z21)
108 *(x23*x23+y23*y23+z23*z23)+1e-6);
109 const REAL angle=(x21*x23+y21*y23+z21*z23)/norm21_norm23;
110 if(angle>=1)
return 0;
111 if(angle<=-1)
return M_PI;
117 const REAL x21=at1.GetX()-at2.GetX();
118 const REAL y21=at1.GetY()-at2.GetY();
119 const REAL z21=at1.GetZ()-at2.GetZ();
121 const REAL x34=at4.GetX()-at3.GetX();
122 const REAL y34=at4.GetY()-at3.GetY();
123 const REAL z34=at4.GetZ()-at3.GetZ();
125 const REAL x23=at3.GetX()-at2.GetX();
126 const REAL y23=at3.GetY()-at2.GetY();
127 const REAL z23=at3.GetZ()-at2.GetZ();
130 const REAL x123= y21*z23-z21*y23;
131 const REAL y123= z21*x23-x21*z23;
132 const REAL z123= x21*y23-y21*x23;
135 const REAL x234= -(y23*z34-z23*y34);
136 const REAL y234= -(z23*x34-x23*z34);
137 const REAL z234= -(x23*y34-y23*x34);
139 const REAL norm123_norm234= sqrt( (x123*x123+y123*y123+z123*z123)
140 *(x234*x234+y234*y234+z234*z234)+1e-6);
142 REAL angle=(x123*x234+y123*y234+z123*z234)/norm123_norm234;
143 if(angle>= 1) angle=0;
146 if(angle<=-1) angle=M_PI;
147 else angle=acos(angle);
149 if((x21*x234 + y21*y234 + z21*z234)<0)
return -angle;
155 const map<
MolAtom*,set<MolAtom*> > &connect,
156 set<MolAtom*> &atomlist,
const MolAtom* finalAtom)
158 const pair<set<MolAtom*>::iterator,
bool> status=atomlist.insert(atom);
159 if(
false==status.second)
return;
160 if(finalAtom==atom)
return;
161 map<MolAtom*,set<MolAtom*> >::const_iterator c=connect.find(atom);
162 set<MolAtom*>::const_iterator pos;
163 for(pos=c->second.begin();pos!=c->second.end();++pos)
170 const map<
MolAtom*,set<MolAtom*> > &connect,
171 map<MolAtom*,unsigned long> &atomlist,
const unsigned long maxdepth,
unsigned long depth)
173 if(atomlist.count(atom)>0)
174 if(atomlist[atom]<=depth)
return;
175 atomlist[atom]=depth;
176 if(depth==maxdepth)
return;
177 map<MolAtom*,set<MolAtom*> >::const_iterator c=connect.find(atom);
178 set<MolAtom*>::const_iterator pos;
179 for(pos=c->second.begin();pos!=c->second.end();++pos)
191 MolAtom::MolAtom(
const REAL x,
const REAL y,
const REAL z,
193 mName(name),mX(x),mY(y),mZ(z),mOccupancy(1.),mpScattPow(pPow),mpMol(&parent),mIsInRing(false),mIsNonFlipAtom(false)
198 VFN_DEBUG_MESSAGE(
"MolAtom::MolAtom()",4)
208 void MolAtom::SetName(
const string &name)
221 catch(
const ObjCrystException &except)
223 cerr<<
"MolAtom::SetName(): Atom parameters not yet declared in a Molecule ?"<<endl;
228 const string& MolAtom::GetName()
const{
return mName;}
229 string& MolAtom::GetName() {
return mName;}
231 const Molecule& MolAtom::GetMolecule()
const{
return *
mpMol;}
232 Molecule& MolAtom::GetMolecule() {
return *
mpMol;}
234 const REAL& MolAtom::X()
const{
return mX;}
235 const REAL& MolAtom::Y()
const{
return mY;}
236 const REAL& MolAtom::Z()
const{
return mZ;}
238 REAL& MolAtom::X(){
return mX;}
239 REAL& MolAtom::Y(){
return mY;}
240 REAL& MolAtom::Z(){
return mZ;}
242 REAL MolAtom::GetX()
const{
return mX;}
243 REAL MolAtom::GetY()
const{
return mY;}
244 REAL MolAtom::GetZ()
const{
return mZ;}
245 REAL MolAtom::GetOccupancy()
const{
return mOccupancy;}
250 void MolAtom::SetOccupancy(
const REAL a){
mOccupancy=a;}
255 void MolAtom::SetScatteringPower(
const ScatteringPower& pow)
264 void MolAtom::XMLOutput(ostream &os,
int indent)
const
266 VFN_DEBUG_ENTRY(
"MolAtom::XMLOutput()",4)
267 for(
int i=0;i<indent;i++) os << " " ;
268 XMLCrystTag tag("Atom",false,true);
269 tag.AddAttribute("Name",this->GetName());
270 if(!this->
IsDummy())tag.AddAttribute("ScattPow",this->GetScatteringPower().GetName());
273 ss.precision(os.precision());
275 tag.AddAttribute(
"x",ss.str());
279 ss.precision(os.precision());
281 tag.AddAttribute(
"y",ss.str());
285 ss.precision(os.precision());
287 tag.AddAttribute(
"z",ss.str());
291 ss.precision(os.precision());
293 tag.AddAttribute(
"Occup",ss.str());
297 VFN_DEBUG_EXIT(
"MolAtom::XMLOutput()",4)
300 void MolAtom::XMLInput(istream &is,
const XMLCrystTag &tag)
302 VFN_DEBUG_ENTRY(
"MolAtom::XMLInput()",7)
304 for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
306 if(
"Name"==tag.GetAttributeName(i))
308 name=tag.GetAttributeValue(i);
310 if(
"ScattPow"==tag.GetAttributeName(i))
314 if(
"x"==tag.GetAttributeName(i))
316 stringstream ss(tag.GetAttributeValue(i));
319 if(
"y"==tag.GetAttributeName(i))
321 stringstream ss(tag.GetAttributeValue(i));
324 if(
"z"==tag.GetAttributeName(i))
326 stringstream ss(tag.GetAttributeValue(i));
329 if(
"Occup"==tag.GetAttributeName(i))
331 stringstream ss(tag.GetAttributeValue(i));
334 if(
"NonFlip"==tag.GetAttributeName(i))
336 stringstream ss(tag.GetAttributeValue(i));
341 VFN_DEBUG_EXIT(
"MolAtom::XMLInput()",7)
345 bool MolAtom::IsInRing()
const{
return mIsInRing;}
364 VFN_DEBUG_ENTRY(
"MolAtom::WXCreate()",5)
366 VFN_DEBUG_EXIT("
MolAtom::WXCreate()",5)
370 void MolAtom::WXDelete(){
if(0!=mpWXCrystObj)
delete mpWXCrystObj;mpWXCrystObj=0;}
371 void MolAtom::WXNotifyDelete(){mpWXCrystObj=0;}
379 const REAL length0,
const REAL sigma,
const REAL delta,
380 Molecule &parent,
const REAL bondOrder):
381 mAtomPair(make_pair(&atom1,&atom2)),
382 mLength0(length0),mDelta(delta),mSigma(sigma),
383 mBondOrder(bondOrder),mIsFreeTorsion(false),
mpMol(&parent)
397 Molecule& MolBond::GetMolecule() {
return *
mpMol;}
400 {
return this->GetAtom1().GetName()+
"-"+this->GetAtom2().GetName();}
402 void MolBond::XMLOutput(ostream &os,
int indent)
const
404 VFN_DEBUG_ENTRY(
"MolBond::XMLOutput()",4)
405 for(
int i=0;i<indent;i++) os << " " ;
407 tag.AddAttribute("Atom1",mAtomPair.first->
GetName());
408 tag.AddAttribute("Atom2",mAtomPair.second->
GetName());
411 ss.precision(os.precision());
413 tag.AddAttribute(
"Length",ss.str());
417 ss.precision(os.precision());
419 tag.AddAttribute(
"Delta",ss.str());
423 ss.precision(os.precision());
425 tag.AddAttribute(
"Sigma",ss.str());
429 ss.precision(os.precision());
431 tag.AddAttribute(
"BondOrder",ss.str());
435 ss.precision(os.precision());
437 tag.AddAttribute(
"FreeTorsion",ss.str());
440 VFN_DEBUG_EXIT(
"MolBond::XMLOutput()",4)
443 void MolBond::XMLInput(istream &is,
const XMLCrystTag &tag)
445 VFN_DEBUG_ENTRY(
"MolBond::XMLInput():",7)
446 for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
448 if(
"Atom1"==tag.GetAttributeName(i))
450 mAtomPair.first=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
452 if(
"Atom2"==tag.GetAttributeName(i))
454 mAtomPair.second=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
456 if(
"Length"==tag.GetAttributeName(i))
458 stringstream ss(tag.GetAttributeValue(i));
461 if(
"Delta"==tag.GetAttributeName(i))
463 stringstream ss(tag.GetAttributeValue(i));
466 if(
"Sigma"==tag.GetAttributeName(i))
468 stringstream ss(tag.GetAttributeValue(i));
471 if(
"BondOrder"==tag.GetAttributeName(i))
473 stringstream ss(tag.GetAttributeValue(i));
476 if(
"FreeTorsion"==tag.GetAttributeName(i))
478 stringstream ss(tag.GetAttributeValue(i));
482 VFN_DEBUG_EXIT(
"MolBond::XMLInput():",7)
489 if(!recalc)
return mLLK;
490 VFN_DEBUG_ENTRY(
"MolBond::GetLogLikelihood():",2)
493 const REAL x=this->GetAtom2().GetX()-this->GetAtom1().GetX();
494 const REAL y=this->GetAtom2().GetY()-this->GetAtom1().GetY();
495 const REAL z=this->GetAtom2().GetZ()-this->GetAtom1().GetZ();
496 const REAL length=sqrt(abs(x*x+y*y+z*z));
500 const REAL tmp2=1/(length+1e-10);
516 mLLK=length-(mLength0+mDelta);
521 #ifdef RESTRAINT_X2_X4_X6
526 mLLK=mLLK2*(1+mLLK2);
531 VFN_DEBUG_EXIT(
"MolBond::GetLogLikelihood():",2)
534 mLLK=length-(mLength0-mDelta);
538 #ifdef RESTRAINT_X2_X4_X6
543 mLLK=mLLK2*(1+mLLK2);
548 VFN_DEBUG_EXIT(
"MolBond::GetLogLikelihood():",2)
561 map<const MolAtom*,XYZ>::const_iterator pos;
562 pos=m.find(mAtomPair.first);
567 pos=m.find(mAtomPair.second);
569 d+= pos->second.x*mDerivAtom2.x
570 +pos->second.y*mDerivAtom2.y
571 +pos->second.z*mDerivAtom2.z;
579 map<MolAtom*,XYZ>::iterator pos;
580 pos=m.find(mAtomPair.first);
587 pos=m.find(mAtomPair.second);
596 cout<<this->
GetName()<<
" :LLK"<<endl;
597 for(map<MolAtom*,XYZ>::const_iterator pos=m.begin();pos!=m.end();++pos)
600 sprintf(buf,
"%10s Grad LLK= (%8.3f %8.3f %8.3f)",
601 pos->first->GetName().c_str(),pos->second.x,pos->second.y,pos->second.z);
607 const MolAtom& MolBond::GetAtom1()
const{
return *(mAtomPair.first);}
608 const MolAtom& MolBond::GetAtom2()
const{
return *(mAtomPair.second);}
609 MolAtom& MolBond::GetAtom1(){
return *(mAtomPair.first);}
610 MolAtom& MolBond::GetAtom2(){
return *(mAtomPair.second);}
611 void MolBond::SetAtom1(MolAtom &at){mAtomPair.first =&at;}
612 void MolBond::SetAtom2(MolAtom &at){mAtomPair.second=&at;}
613 REAL MolBond::GetLength()
const
618 REAL MolBond::GetLength0()
const{
return mLength0;}
619 REAL MolBond::GetLengthDelta()
const{
return mDelta;}
620 REAL MolBond::GetLengthSigma()
const{
return mSigma;}
621 REAL MolBond::GetBondOrder()
const{
return mBondOrder;}
623 REAL& MolBond::Length0(){
return mLength0;}
624 REAL& MolBond::LengthDelta(){
return mDelta;}
625 REAL& MolBond::LengthSigma(){
return mSigma;}
626 REAL& MolBond::BondOrder(){
return mBondOrder;}
628 void MolBond::SetLength0(
const REAL a){mLength0=a;}
629 void MolBond::SetLengthDelta(
const REAL a){mDelta=a;}
630 void MolBond::SetLengthSigma(
const REAL a){mSigma=a;}
631 void MolBond::SetBondOrder(
const REAL a){mBondOrder=a;}
633 bool MolBond::IsFreeTorsion()
const{
return mIsFreeTorsion;}
634 void MolBond::SetFreeTorsion(
const bool isFreeTorsion)
636 if(mIsFreeTorsion==isFreeTorsion)
return;
637 mIsFreeTorsion=isFreeTorsion;
646 VFN_DEBUG_ENTRY(
"MolBond::WXCreate()",5)
648 VFN_DEBUG_EXIT("
MolBond::WXCreate()",5)
652 void MolBond::WXDelete(){
if(0!=mpWXCrystObj)
delete mpWXCrystObj;mpWXCrystObj=0;}
653 void MolBond::WXNotifyDelete(){mpWXCrystObj=0;}
661 const REAL angle,
const REAL sigma,
const REAL delta,
663 mAngle0(angle),mDelta(delta),mSigma(sigma),
mpMol(&parent)
680 const Molecule& MolBondAngle::GetMolecule()
const{
return *
mpMol;}
681 Molecule& MolBondAngle::GetMolecule() {
return *
mpMol;}
683 string MolBondAngle::GetName()
const
685 return this->GetAtom1().GetName()+
"-"
686 +this->GetAtom2().GetName()+
"-"
687 +this->GetAtom3().GetName();
690 void MolBondAngle::XMLOutput(ostream &os,
int indent)
const
692 VFN_DEBUG_ENTRY(
"MolBondAngle::XMLOutput()",4)
693 for(
int i=0;i<indent;i++) os << " " ;
694 XMLCrystTag tag("BondAngle",false,true);
695 tag.AddAttribute("Atom1",this->GetAtom1().GetName());
696 tag.AddAttribute("Atom2",this->GetAtom2().GetName());
697 tag.AddAttribute("Atom3",this->GetAtom3().GetName());
700 ss.precision(os.precision());
701 ss <<mAngle0*RAD2DEG;
702 tag.AddAttribute(
"Angle",ss.str());
706 ss.precision(os.precision());
708 tag.AddAttribute(
"Delta",ss.str());
712 ss.precision(os.precision());
714 tag.AddAttribute(
"Sigma",ss.str());
717 VFN_DEBUG_EXIT(
"MolBondAngle::XMLOutput()",4)
720 void MolBondAngle::XMLInput(istream &is,
const XMLCrystTag &tag)
722 VFN_DEBUG_ENTRY(
"MolBondAngle::XMLInput():",4)
724 for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
726 if(
"Atom1"==tag.GetAttributeName(i))
730 if(
"Atom2"==tag.GetAttributeName(i))
734 if(
"Atom3"==tag.GetAttributeName(i))
738 if(
"Angle"==tag.GetAttributeName(i))
740 stringstream ss(tag.GetAttributeValue(i));
744 if(
"Delta"==tag.GetAttributeName(i))
746 stringstream ss(tag.GetAttributeValue(i));
750 if(
"Sigma"==tag.GetAttributeName(i))
752 stringstream ss(tag.GetAttributeValue(i));
757 VFN_DEBUG_EXIT(
"MolBondAngle::XMLInput():",4)
759 REAL& MolBondAngle::Angle0()
763 REAL& MolBondAngle::AngleDelta(){
return mDelta;}
764 REAL& MolBondAngle::AngleSigma(){
return mSigma;}
766 REAL MolBondAngle::GetAngle0()
const{
return mAngle0;}
767 REAL MolBondAngle::GetAngleDelta()
const{
return mDelta;}
768 REAL MolBondAngle::GetAngleSigma()
const{
return mSigma;}
770 void MolBondAngle::SetAngle0(
const REAL angle){mAngle0=angle;}
771 void MolBondAngle::SetAngleDelta(
const REAL delta){mDelta=delta;}
772 void MolBondAngle::SetAngleSigma(
const REAL sigma){mSigma=sigma;}
774 REAL MolBondAngle::GetAngle()
const
776 return GetBondAngle(this->GetAtom1(),this->GetAtom2(),this->GetAtom3());
783 if(!recalc)
return mLLK;
784 VFN_DEBUG_ENTRY(
"MolBondAngle::GetLogLikelihood():",2)
787 const REAL x21=this->GetAtom1().GetX()-this->GetAtom2().GetX();
788 const REAL y21=this->GetAtom1().GetY()-this->GetAtom2().GetY();
789 const REAL z21=this->GetAtom1().GetZ()-this->GetAtom2().GetZ();
790 const REAL x23=this->GetAtom3().GetX()-this->GetAtom2().GetX();
791 const REAL y23=this->GetAtom3().GetY()-this->GetAtom2().GetY();
792 const REAL z23=this->GetAtom3().GetZ()-this->GetAtom2().GetZ();
794 const REAL n1=sqrt(abs(x21*x21+y21*y21+z21*z21));
795 const REAL n3=sqrt(abs(x23*x23+y23*y23+z23*z23));
796 const REAL p=x21*x23+y21*y23+z21*z23;
798 const REAL a0=p/(n1*n3+1e-10);
803 if(a0<=-1) angle=M_PI;
809 const REAL s=1/(sqrt(1-a0*a0+1e-6));
810 const REAL s0=-s/(n1*n3+1e-10);
811 const REAL s1= s*p/(n3*n1*n1*n1+1e-10);
812 const REAL s3= s*p/(n1*n3*n3*n3+1e-10);
817 mDerivAtom3.x=s0*x21+s3*x23;
818 mDerivAtom3.y=s0*y21+s3*y23;
819 mDerivAtom3.z=s0*z21+s3*z23;
833 mLLK=angle-(mAngle0+mDelta);
837 #ifdef RESTRAINT_X2_X4_X6
842 mLLK=mLLK2*(1+mLLK2);
847 VFN_DEBUG_EXIT(
"MolBondAngle::GetLogLikelihood():",2)
850 mLLK=angle-(mAngle0-mDelta);
854 #ifdef RESTRAINT_X2_X4_X6
859 mLLK=mLLK2*(1+mLLK2);
864 VFN_DEBUG_EXIT(
"MolBondAngle::GetLogLikelihood():",2)
877 map<const MolAtom*,XYZ>::const_iterator pos;
885 d+= pos->second.x*mDerivAtom2.x
886 +pos->second.y*mDerivAtom2.y
887 +pos->second.z*mDerivAtom2.z;
890 d+= pos->second.x*mDerivAtom3.x
891 +pos->second.y*mDerivAtom3.y
892 +pos->second.z*mDerivAtom3.z;
900 map<MolAtom*,XYZ>::iterator pos;
924 cout<<this->GetName()<<
" :LLK"<<endl;
925 for(map<MolAtom*,XYZ>::const_iterator pos=m.begin();pos!=m.end();++pos)
928 sprintf(buf,
"%10s Grad LLK= (%8.3f %8.3f %8.3f)",
929 pos->first->GetName().c_str(),pos->second.x,pos->second.y,pos->second.z);
935 const MolAtom& MolBondAngle::GetAtom1()
const{
return *(
mvpAtom[0]);}
936 const MolAtom& MolBondAngle::GetAtom2()
const{
return *(
mvpAtom[1]);}
937 const MolAtom& MolBondAngle::GetAtom3()
const{
return *(
mvpAtom[2]);}
938 MolAtom& MolBondAngle::GetAtom1(){
return *(
mvpAtom[0]);}
939 MolAtom& MolBondAngle::GetAtom2(){
return *(
mvpAtom[1]);}
940 MolAtom& MolBondAngle::GetAtom3(){
return *(
mvpAtom[2]);}
941 void MolBondAngle::SetAtom1(MolAtom& at){
mvpAtom[0]=&at;}
942 void MolBondAngle::SetAtom2(MolAtom& at){
mvpAtom[1]=&at;}
943 void MolBondAngle::SetAtom3(MolAtom& at){
mvpAtom[2]=&at;}
947 std::size_t MolBondAngle::size()
const {
return mvpAtom.size();}
948 vector<MolAtom*>::const_iterator MolBondAngle::begin()
const {
return mvpAtom.begin();}
949 vector<MolAtom*>::const_iterator MolBondAngle::end()
const {
return mvpAtom.end();}
956 VFN_DEBUG_ENTRY(
"MolBondAngle::WXCreate()",5)
962 void MolBondAngle::WXDelete(){
if(0!=mpWXCrystObj)
delete mpWXCrystObj;mpWXCrystObj=0;}
963 void MolBondAngle::WXNotifyDelete(){mpWXCrystObj=0;}
972 const REAL angle,
const REAL sigma,
const REAL delta,
974 mAngle0(angle),mDelta(delta),mSigma(sigma),
mpMol(&parent)
979 VFN_DEBUG_ENTRY(
"MolDihedralAngle::MolDihedralAngle()",5)
985 mAngle0=fmod((REAL)mAngle0,(REAL)(2*M_PI));
986 if(mAngle0<(-M_PI)) mAngle0+=2*M_PI;
987 if(mAngle0>M_PI) mAngle0-=2*M_PI;
988 VFN_DEBUG_EXIT(
"MolDihedralAngle::MolDihedralAngle()",5)
998 const Molecule& MolDihedralAngle::GetMolecule()
const{
return *
mpMol;}
999 Molecule& MolDihedralAngle::GetMolecule() {
return *
mpMol;}
1001 string MolDihedralAngle::GetName()
const
1003 return this->GetAtom1().GetName()+
"-"
1004 +this->GetAtom2().GetName()+
"-"
1005 +this->GetAtom3().GetName()+
"-"
1006 +this->GetAtom4().GetName();
1009 void MolDihedralAngle::XMLOutput(ostream &os,
int indent)
const
1011 VFN_DEBUG_ENTRY(
"MolDihedralAngle::XMLOutput()",4)
1012 for(
int i=0;i<indent;i++) os << " " ;
1013 XMLCrystTag tag("DihedralAngle",false,true);
1014 tag.AddAttribute("Atom1",this->GetAtom1().GetName());
1015 tag.AddAttribute("Atom2",this->GetAtom2().GetName());
1016 tag.AddAttribute("Atom3",this->GetAtom3().GetName());
1017 tag.AddAttribute("Atom4",this->GetAtom4().GetName());
1020 ss.precision(os.precision());
1021 ss <<mAngle0*RAD2DEG;
1022 tag.AddAttribute(
"Angle",ss.str());
1026 ss.precision(os.precision());
1027 ss <<mDelta*RAD2DEG;
1028 tag.AddAttribute(
"Delta",ss.str());
1032 ss.precision(os.precision());
1033 ss <<mSigma*RAD2DEG;
1034 tag.AddAttribute(
"Sigma",ss.str());
1037 VFN_DEBUG_EXIT(
"MolDihedralAngle::XMLOutput()",4)
1040 void MolDihedralAngle::XMLInput(istream &is,
const XMLCrystTag &tag)
1042 VFN_DEBUG_ENTRY(
"MolDihedralAngle::XMLInput():",5)
1044 for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
1046 if(
"Atom1"==tag.GetAttributeName(i))
1048 mvpAtom[0]=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
1050 if(
"Atom2"==tag.GetAttributeName(i))
1052 mvpAtom[1]=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
1054 if(
"Atom3"==tag.GetAttributeName(i))
1056 mvpAtom[2]=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
1058 if(
"Atom4"==tag.GetAttributeName(i))
1060 mvpAtom[3]=&(
mpMol->GetAtom(tag.GetAttributeValue(i)));
1062 if(
"Angle"==tag.GetAttributeName(i))
1064 stringstream ss(tag.GetAttributeValue(i));
1068 if(
"Delta"==tag.GetAttributeName(i))
1070 stringstream ss(tag.GetAttributeValue(i));
1074 if(
"Sigma"==tag.GetAttributeName(i))
1076 stringstream ss(tag.GetAttributeValue(i));
1081 VFN_DEBUG_EXIT(
"MolDihedralAngle::XMLInput():",5)
1084 REAL MolDihedralAngle::GetAngle()
const
1087 const REAL angle=
GetDihedralAngle(this->GetAtom1(),this->GetAtom2(),this->GetAtom3(),this->GetAtom4());
1088 if((angle-mAngle0)>M_PI)
return angle-2*M_PI;
1089 else if((angle-mAngle0)<(-M_PI))
return angle+2*M_PI;
1094 REAL& MolDihedralAngle::Angle0(){
return mAngle0;}
1095 REAL& MolDihedralAngle::AngleDelta(){
return mDelta;}
1096 REAL& MolDihedralAngle::AngleSigma(){
return mSigma;}
1098 REAL MolDihedralAngle::GetAngle0()
const{
return mAngle0;}
1099 REAL MolDihedralAngle::GetAngleDelta()
const{
return mDelta;}
1100 REAL MolDihedralAngle::GetAngleSigma()
const{
return mSigma;}
1102 void MolDihedralAngle::SetAngle0(
const REAL angle){mAngle0=angle;}
1103 void MolDihedralAngle::SetAngleDelta(
const REAL delta){mDelta=delta;}
1104 void MolDihedralAngle::SetAngleSigma(
const REAL sigma){mSigma=sigma;}
1110 if(!recalc)
return mLLK;
1111 VFN_DEBUG_ENTRY(
"MolDihedralAngle::GetLogLikelihood():",2)
1113 const REAL x21=this->GetAtom1().GetX()-this->GetAtom2().GetX();
1114 const REAL y21=this->GetAtom1().GetY()-this->GetAtom2().GetY();
1115 const REAL z21=this->GetAtom1().GetZ()-this->GetAtom2().GetZ();
1117 const REAL x34=this->GetAtom4().GetX()-this->GetAtom3().GetX();
1118 const REAL y34=this->GetAtom4().GetY()-this->GetAtom3().GetY();
1119 const REAL z34=this->GetAtom4().GetZ()-this->GetAtom3().GetZ();
1121 const REAL x23=this->GetAtom3().GetX()-this->GetAtom2().GetX();
1122 const REAL y23=this->GetAtom3().GetY()-this->GetAtom2().GetY();
1123 const REAL z23=this->GetAtom3().GetZ()-this->GetAtom2().GetZ();
1126 const REAL x123= y21*z23-z21*y23;
1127 const REAL y123= z21*x23-x21*z23;
1128 const REAL z123= x21*y23-y21*x23;
1131 const REAL x234= -(y23*z34-z23*y34);
1132 const REAL y234= -(z23*x34-x23*z34);
1133 const REAL z234= -(x23*y34-y23*x34);
1135 const REAL n123= sqrt(x123*x123+y123*y123+z123*z123+1e-7);
1136 const REAL n234= sqrt(x234*x234+y234*y234+z234*z234+1e-6);
1138 const REAL p=x123*x234+y123*y234+z123*z234;
1139 const REAL a0=p/(n123*n234+1e-10);
1144 if(a0<=-1) angle=M_PI;
1145 else angle=acos(a0);
1148 if((x21*x34 + y21*y34 + z21*z34)<0) {angle=-angle;sgn=-1;}
1153 const REAL s=sgn/(sqrt(1-a0*a0+1e-6));
1154 const REAL s0=-s/(n123*n234+1e-10);
1155 const REAL s1= s*p/(n234*n123*n123*n123+1e-10);
1156 const REAL s3= s*p/(n123*n234*n234*n234+1e-10);
1157 mDerivAtom1.x=s0*(-z23*y234+y23*z234)+s1*(-z23*y123+y23*z123);
1158 mDerivAtom1.y=s0*(-x23*z234+z23*x234)+s1*(-x23*z123+z23*x123);
1159 mDerivAtom1.z=s0*(-y23*x234+x23*y234)+s1*(-y23*x123+x23*y123);
1161 mDerivAtom4.x=s0*(-z23*y123+y23*z123)+s3*(-z23*y234+y23*z234);
1162 mDerivAtom4.y=s0*(-x23*z123+z23*x123)+s3*(-x23*z234+z23*x234);
1163 mDerivAtom4.z=s0*(-y23*x123+x23*y123)+s3*(-y23*x234+x23*y234);
1165 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);
1166 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);
1167 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);
1169 mDerivAtom3.x=-(
mDerivAtom1.x+mDerivAtom2.x+mDerivAtom4.x);
1170 mDerivAtom3.y=-(
mDerivAtom1.y+mDerivAtom2.y+mDerivAtom4.y);
1171 mDerivAtom3.z=-(
mDerivAtom1.z+mDerivAtom2.z+mDerivAtom4.z);
1180 mLLK=angle-(mAngle0+mDelta);
1186 #ifdef RESTRAINT_X2_X4_X6
1191 mLLK=mLLK2*(1+mLLK2);
1196 VFN_DEBUG_EXIT(
"MolDihedralAngle::GetLogLikelihood():",2)
1199 mLLK=angle-(mAngle0-mDelta);
1205 #ifdef RESTRAINT_X2_X4_X6
1210 mLLK=mLLK2*(1+mLLK2);
1215 VFN_DEBUG_EXIT(
"MolDihedralAngle::GetLogLikelihood():",2)
1228 map<const MolAtom*,XYZ>::const_iterator pos;
1236 d+= pos->second.x*mDerivAtom2.x
1237 +pos->second.y*mDerivAtom2.y
1238 +pos->second.z*mDerivAtom2.z;
1241 d+= pos->second.x*mDerivAtom3.x
1242 +pos->second.y*mDerivAtom3.y
1243 +pos->second.z*mDerivAtom3.z;
1246 d+= pos->second.x*mDerivAtom4.x
1247 +pos->second.y*mDerivAtom4.y
1248 +pos->second.z*mDerivAtom4.z;
1256 map<MolAtom*,XYZ>::iterator pos;
1287 cout<<this->GetName()<<
" :LLK"<<endl;
1288 for(map<MolAtom*,XYZ>::const_iterator pos=m.begin();pos!=m.end();++pos)
1291 sprintf(buf,
"%10s Grad LLK= (%8.3f %8.3f %8.3f)",
1292 pos->first->GetName().c_str(),pos->second.x,pos->second.y,pos->second.z);
1298 const MolAtom& MolDihedralAngle::GetAtom1()
const{
return *(
mvpAtom[0]);}
1299 const MolAtom& MolDihedralAngle::GetAtom2()
const{
return *(
mvpAtom[1]);}
1300 const MolAtom& MolDihedralAngle::GetAtom3()
const{
return *(
mvpAtom[2]);}
1301 const MolAtom& MolDihedralAngle::GetAtom4()
const{
return *(
mvpAtom[3]);}
1302 void MolDihedralAngle::SetAtom1(MolAtom& at){
mvpAtom[0]=&at;}
1303 void MolDihedralAngle::SetAtom2(MolAtom& at){
mvpAtom[1]=&at;}
1304 void MolDihedralAngle::SetAtom3(MolAtom& at){
mvpAtom[2]=&at;}
1305 void MolDihedralAngle::SetAtom4(MolAtom& at){
mvpAtom[3]=&at;}
1306 MolAtom& MolDihedralAngle::GetAtom1(){
return *(
mvpAtom[0]);}
1307 MolAtom& MolDihedralAngle::GetAtom2(){
return *(
mvpAtom[1]);}
1308 MolAtom& MolDihedralAngle::GetAtom3(){
return *(
mvpAtom[2]);}
1309 MolAtom& MolDihedralAngle::GetAtom4(){
return *(
mvpAtom[3]);}
1310 std::size_t MolDihedralAngle::size()
const {
return mvpAtom.size();}
1311 vector<MolAtom*>::const_iterator MolDihedralAngle::begin()
const {
return mvpAtom.begin();}
1312 vector<MolAtom*>::const_iterator MolDihedralAngle::end()
const {
return mvpAtom.end();}
1316 #ifdef __WX__CRYST__
1319 VFN_DEBUG_ENTRY(
"MolDihedralAngle::WXCreate()",5)
1322 return mpWXCrystObj;
1325 void MolDihedralAngle::WXDelete(){
if(0!=mpWXCrystObj)
delete mpWXCrystObj;mpWXCrystObj=0;}
1326 void MolDihedralAngle::WXNotifyDelete(){mpWXCrystObj=0;}
1333 string RigidGroup::GetName()
const
1335 set<MolAtom *>::const_iterator at=this->begin();
1336 string name=(*at++)->GetName();
1337 for(;at!=this->end();++at) name+=
", "+(*at)->GetName();
1351 const std::list<MolAtom*>& MolRing::GetAtomList()
const
1354 std::list<MolAtom*>& MolRing::GetAtomList()
1365 mQ0(1),mQ1(0),mQ2(0),mQ3(0),mIsUniQuaternion(true)
1367 VFN_DEBUG_MESSAGE(
"Quaternion::Quaternion()",5)
1375 mQ0(q0),mQ1(q1),mQ2(q2),mQ3(q3),mIsUniQuaternion(unit)
1377 VFN_DEBUG_MESSAGE(
"Quaternion::Quaternion()",5)
1381 Quaternion::~Quaternion()
1383 VFN_DEBUG_MESSAGE(
"Quaternion::~Quaternion()",5)
1391 VFN_DEBUG_MESSAGE(
"Quaternion::RotationQuaternion()",4)
1392 const REAL s=sin(ang/2.)/sqrt(v1*v1+v2*v2+v3*v3+1e-7);
1393 return Quaternion(cos(ang/2.),s*v1,s*v2,s*v3,
1405 (this->Q0()*q.Q0()-this->Q1()*q.Q1()-this->Q2()*q.Q2()-this->Q3()*q.Q3(),
1406 this->Q0()*q.Q1()+this->Q1()*q.Q0()+this->Q2()*q.Q3()-this->Q3()*q.Q2(),
1407 this->Q0()*q.Q2()-this->Q1()*q.Q3()+this->Q2()*q.Q0()+this->Q3()*q.Q1(),
1408 this->Q0()*q.Q3()+this->Q1()*q.Q2()-this->Q2()*q.Q1()+this->Q3()*q.Q0(),
false);
1411 void Quaternion::operator*=(
const Quaternion &q)
1415 const REAL q1=this->Q0()*q.Q1()+this->Q1()*q.Q0()+this->Q2()*q.Q3()-this->Q3()*q.Q2();
1416 const REAL q2=this->Q0()*q.Q2()+this->Q2()*q.Q0()-this->Q1()*q.Q3()+this->Q3()*q.Q1();
1417 const REAL q3=this->Q0()*q.Q3()+this->Q3()*q.Q0()+this->Q1()*q.Q2()-this->Q2()*q.Q1();
1418 this->Q0()= this->Q0()*q.Q0()-this->Q1()*q.Q1()-this->Q2()*q.Q2()-this->Q3()*q.Q3();
1426 void Quaternion::XMLOutput(ostream &os,
int indent)
const
1428 VFN_DEBUG_ENTRY(
"Quaternion::XMLOutput()",4)
1429 for(
int i=0;i<indent;i++) os << " " ;
1434 ss.precision(os.precision());
1436 tag.AddAttribute(
"Q0",ss.str());
1440 ss.precision(os.precision());
1442 tag.AddAttribute(
"Q1",ss.str());
1446 ss.precision(os.precision());
1448 tag.AddAttribute(
"Q2",ss.str());
1452 ss.precision(os.precision());
1454 tag.AddAttribute(
"Q3",ss.str());
1458 ss.precision(os.precision());
1459 ss <<mIsUniQuaternion;
1460 tag.AddAttribute(
"IsUnitQuaternion",ss.str());
1463 VFN_DEBUG_EXIT(
"Quaternion::XMLOutput()",4)
1466 void Quaternion::XMLInput(istream &is,
const XMLCrystTag &tag)
1468 VFN_DEBUG_ENTRY(
"Quaternion::XMLInput()",5)
1469 for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
1471 if(
"Q0"==tag.GetAttributeName(i))
1473 stringstream ss(tag.GetAttributeValue(i));
1476 if(
"Q1"==tag.GetAttributeName(i))
1478 stringstream ss(tag.GetAttributeValue(i));
1481 if(
"Q2"==tag.GetAttributeName(i))
1483 stringstream ss(tag.GetAttributeValue(i));
1486 if(
"Q3"==tag.GetAttributeName(i))
1488 stringstream ss(tag.GetAttributeValue(i));
1491 if(
"IsUnitQuaternion"==tag.GetAttributeName(i))
1493 stringstream ss(tag.GetAttributeValue(i));
1494 ss >>mIsUniQuaternion;
1498 VFN_DEBUG_EXIT(
"Quaternion::XMLInput()",5)
1514 const REAL p0=-mQ1*v1 - mQ2*v2 - mQ3*v3;
1515 const REAL p1=
mQ0*v1 + mQ2*v3 - mQ3*v2;
1516 const REAL p2=
mQ0*v2 - mQ1*v3 + mQ3*v1;
1517 const REAL p3=
mQ0*v3 + mQ1*v2 - mQ2*v1;
1519 v1 =-p0*mQ1 + p1*
mQ0 - p2*mQ3 + p3*mQ2;
1520 v2 =-p0*mQ2 + p2*
mQ0 + p1*mQ3 - p3*mQ1;
1521 v3 =-p0*mQ3 + p3*
mQ0 - p1*mQ2 + p2*mQ1;
1526 const REAL norm=sqrt( this->Q0()*this->Q0()
1527 +this->Q1()*this->Q1()
1528 +this->Q2()*this->Q2()
1529 +this->Q3()*this->Q3());
1535 REAL Quaternion::GetNorm()
const
1536 {
return sqrt( this->Q0()*this->Q0()
1537 +this->Q1()*this->Q1()
1538 +this->Q2()*this->Q2()
1539 +this->Q3()*this->Q3());}
1541 const REAL& Quaternion::Q0()
const{
return mQ0;}
1542 const REAL& Quaternion::Q1()
const{
return mQ1;}
1543 const REAL& Quaternion::Q2()
const{
return mQ2;}
1544 const REAL& Quaternion::Q3()
const{
return mQ3;}
1545 REAL& Quaternion::Q0(){
return mQ0;}
1546 REAL& Quaternion::Q1(){
return mQ1;}
1547 REAL& Quaternion::Q2(){
return mQ2;}
1548 REAL& Quaternion::Q3(){
return mQ3;}
1554 StretchMode::~StretchMode(){}
1558 mpAtom0(&at0),mpAtom1(&at1),mpBond(pBond)
1561 mpMol = &(at1.GetMolecule());
1564 StretchModeBondLength::~StretchModeBondLength(){}
1574 const REAL n=sqrt(dx*dx+dy*dy+dz*dz+1e-7);
1591 for(map<const MolBond*,REAL>::const_iterator pos=this->
mvpBrokenBond.begin();
1593 for(map<const MolBondAngle*,REAL>::const_iterator pos=this->
mvpBrokenBondAngle.begin();
1605 cout<<
", Translated Atoms:";
1609 cout<<(*atom)->GetName()<<
",";
1615 const bool keepCenter)
1620 const REAL l=sqrt(dx*dx+dy*dy+dz*dz+1e-7);
1622 const REAL change=amplitude/l;
1630 const bool keepCenter)
1637 mpAtom0(&at0),mpAtom1(&at1),mpAtom2(&at2),mpBondAngle(pBondAngle)
1640 mpMol = &(at1.GetMolecule());
1643 StretchModeBondAngle::~StretchModeBondAngle(){}
1649 const REAL x1=
mpAtom1->GetX(),
1653 const REAL dx10=
mpAtom0->GetX()-x1,
1660 REAL vx=dy10*dz12-dz10*dy12,
1661 vy=dz10*dx12-dx10*dz12,
1662 vz=dx10*dy12-dy10*dx12;
1664 const REAL n=sqrt(vx*vx+vy*vy+vz*vz+1e-10);
1679 const REAL x=(*pos)->GetX()-x1,
1680 y=(*pos)->GetY()-y1,
1681 z=(*pos)->GetZ()-z1;
1690 for(map<const MolBond*,REAL>::const_iterator pos=this->
mvpBrokenBond.begin();
1692 for(map<const MolBondAngle*,REAL>::const_iterator pos=this->
mvpBrokenBondAngle.begin();
1704 cout<<
", Rotated Atoms:";
1708 cout<<(*atom)->GetName()<<
",";
1714 const bool keepCenter)
1723 const REAL vx=dy10*dz12-dz10*dy12;
1724 const REAL vy=dz10*dx12-dx10*dz12;
1725 const REAL vz=dx10*dy12-dy10*dx12;
1730 const bool keepCenter)
1737 mpAtom1(&at1),mpAtom2(&at2),mpDihedralAngle(pAngle)
1740 mpMol = &(at1.GetMolecule());
1743 StretchModeTorsion::~StretchModeTorsion(){}
1750 const REAL x1=
mpAtom1->GetX(),
1758 const REAL n=sqrt(vx*vx+vy*vy+vz*vz+1e-10);
1767 const REAL x=(*pos)->GetX()-x1,
1768 y=(*pos)->GetY()-y1,
1769 z=(*pos)->GetZ()-z1;
1778 for(map<const MolBond*,REAL>::const_iterator pos=this->
mvpBrokenBond.begin();
1780 for(map<const MolBondAngle*,REAL>::const_iterator pos=this->
mvpBrokenBondAngle.begin();
1792 cout<<
", Rotated Atoms:";
1796 cout<<(*atom)->GetName()<<
",";
1807 const bool keepCenter)
1819 mpAtom1(&at1),mpAtom2(&at2)
1822 mpMol = &(at1.GetMolecule());
1825 StretchModeTwist::~StretchModeTwist(){}
1831 const REAL x1=
mpAtom1->GetX(),
1839 const REAL n=sqrt(vx*vx+vy*vy+vz*vz+1e-10);
1848 const REAL x=(*pos)->GetX()-x1,
1849 y=(*pos)->GetY()-y1,
1850 z=(*pos)->GetZ()-z1;
1859 for(map<const MolBond*,REAL>::const_iterator pos=this->
mvpBrokenBond.begin();
1861 for(map<const MolBondAngle*,REAL>::const_iterator pos=this->
mvpBrokenBondAngle.begin();
1873 os<<
", Rotated Atoms:";
1877 os<<(*atom)->GetName()<<
",";
1888 const bool keepCenter)
1893 if((abs(dx)+abs(dy)+abs(dz))<1e-6)
return;
1894 const REAL change=(REAL)(2.*rand()-RAND_MAX)/(REAL)RAND_MAX*
mBaseAmplitude*amplitude;
1907 set<MolBondAngle*> &va,
1908 set<MolDihedralAngle*> &vd):
1913 for(set<MolBond*>::iterator pos=vb.begin();pos!=vb.end();++pos)
1914 mvpBond.push_back(*pos);
1915 for(set<MolBondAngle*>::iterator pos=va.begin();pos!=va.end();++pos)
1916 mvpBondAngle.push_back(*pos);
1917 for(set<MolDihedralAngle*>::iterator pos=vd.begin();pos!=vd.end();++pos)
1918 mvpDihedralAngle.push_back(*pos);
1923 if(full) os<<
"MDAtomGroup: ";
1924 for(set<MolAtom*>::const_iterator pos=mvpAtom.begin();pos!=mvpAtom.end();++pos)
1925 os<<(*pos)->GetName()<<
" ";
1928 os<<endl<<
" Involving bond restraints:";
1929 for(vector<MolBond*>::const_iterator pos=mvpBond.begin();pos!=mvpBond.end();++pos)
1930 os<<(*pos)->GetName()<<
" ";
1931 os<<endl<<
" Involving bond angle restraints:";
1932 for(vector<MolBondAngle*>::const_iterator pos=mvpBondAngle.begin();pos!=mvpBondAngle.end();++pos)
1933 os<<(*pos)->GetName()<<
" ";
1934 os<<endl<<
" Involving dihedral angle restraints:";
1935 for(vector<MolDihedralAngle*>::const_iterator pos=mvpDihedralAngle.begin();pos!=mvpDihedralAngle.end();++pos)
1936 os<<(*pos)->GetName()<<
" ";
1947 mDeleteSubObjInDestructor(1), mBaseRotationAmplitude(M_PI*0.02), mIsSelfOptimizing(false),
1948 mpCenterAtom(0), mMDMoveFreq(0.0), mMDMoveEnergy(40.), mLogLikelihoodScale(1.0)
1950 VFN_DEBUG_MESSAGE(
"Molecule::Molecule()",5)
1954 gpRefParTypeScattTranslX,
1955 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
1962 gpRefParTypeScattTranslY,
1963 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
1970 gpRefParTypeScattTranslZ,
1971 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
1978 gpRefParTypeScattOccup,
1979 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.,1.);
1986 gpRefParTypeScattOrient,
1987 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
1995 gpRefParTypeScattOrient,
1996 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
2004 gpRefParTypeScattOrient,
2005 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
2013 gpRefParTypeScattOrient,
2014 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
2021 mLocalParamSet=this->
CreateParamSet(
"saved parameters for local minimization");
2035 mDeleteSubObjInDestructor(old.mDeleteSubObjInDestructor), mIsSelfOptimizing(false), mpCenterAtom(0)
2037 VFN_DEBUG_ENTRY(
"Molecule::Molecule(old&)",5)
2042 gpRefParTypeScattTranslX,
2043 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
2050 gpRefParTypeScattTranslY,
2051 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
2058 gpRefParTypeScattTranslZ,
2059 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
true,1.,1.);
2066 gpRefParTypeScattOccup,
2067 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.,1.);
2074 gpRefParTypeScattOrient,
2075 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
2083 gpRefParTypeScattOrient,
2084 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
2092 gpRefParTypeScattOrient,
2093 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
2101 gpRefParTypeScattOrient,
2102 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
2108 mLocalParamSet=this->
CreateParamSet(
"saved parameters for local minimization");
2131 VFN_DEBUG_EXIT(
"Molecule::Molecule(old&)",5)
2136 VFN_DEBUG_ENTRY(
"Molecule::~Molecule()",5)
2140 vector<MolAtom*>::iterator pos;
2144 vector<MolBond*>::iterator pos;
2148 vector<MolBondAngle*>::iterator pos;
2152 vector<MolDihedralAngle*>::iterator pos;
2156 VFN_DEBUG_EXIT(
"Molecule::~Molecule()",5)
2161 VFN_DEBUG_MESSAGE(
"Molecule::CreateCopy()",5)
2167 static const string className=
"Molecule";
2173 if(
mName==name)
return;
2189 cerr<<
"Molecule::SetName(): parameters not yet declared in a Molecule ?"<<endl;
2196 std::map<std::string,float> velts;
2197 for(std::vector<MolAtom*>::const_iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
2199 if((*pos)->IsDummy())
continue;
2200 string p=(*pos)->GetScatteringPower().GetSymbol();
2201 if(velts.count(p)==0)
2202 velts[p]=(*pos)->GetOccupancy();
2203 else velts[p]+=(*pos)->GetOccupancy();
2206 s<<std::setprecision(2);
2207 for(std::map<std::string,float>::const_iterator pos=velts.begin();pos!=velts.end();++pos)
2209 if(pos!=velts.begin()) s<<
" ";
2210 float nb=pos->second;
2211 if(abs(round(nb)-nb)<0.005)
2213 if(
int(round(nb))==1) s<<pos->first;
2214 else s<<pos->first<<int(round(nb));
2216 else s<<pos->first<<nb;
2223 VFN_DEBUG_MESSAGE(
"Molecule::Print()",5)
2229 VFN_DEBUG_ENTRY(
"Molecule::XMLOutput()",4)
2234 for(
int i=0;i<indent;i++) os <<
" " ;
2236 tag.AddAttribute(
"Name",
mName);
2239 tag.AddAttribute(
"MDMoveFreq",sst.str());
2242 tag.AddAttribute(
"MDMoveEnergy",sst.str());
2245 tag.AddAttribute(
"LogLikelihoodScale",sst.str());
2250 mQuat.XMLOutput(os,indent);
2272 vector<MolAtom*>::const_iterator pos;
2274 (*pos)->XMLOutput(os,indent);
2277 vector<MolBond*>::const_iterator pos;
2279 (*pos)->XMLOutput(os,indent);
2282 vector<MolBondAngle*>::const_iterator pos;
2284 (*pos)->XMLOutput(os,indent);
2287 vector<MolDihedralAngle*>::const_iterator pos;
2289 (*pos)->XMLOutput(os,indent);
2292 vector<RigidGroup*>::const_iterator pos;
2300 for (set<MolAtom*>::const_iterator at = (*pos)->begin(); at != (*pos)->end(); ++at)
2301 tagg.AddAttribute((boost::format(
"Atom%d") %idx++).str(), (*at)->GetName());
2311 for(
int i=0;i<indent;i++) os <<
" " ;
2320 for(
int i=0;i<indent;i++) os <<
" " ;
2325 tag.SetIsEndTag(
true);
2326 for(
int i=0;i<indent;i++) os <<
" " ;
2328 VFN_DEBUG_EXIT(
"Molecule::XMLOutput()",4)
2333 VFN_DEBUG_ENTRY(
"Molecule::XMLInput()",5)
2334 for(
unsigned int i=0;i<tag.GetNbAttribute();i++)
2336 if(
"Name"==tag.GetAttributeName(i))
2338 mName=tag.GetAttributeValue(i);
2340 if(
"MDMoveFreq"==tag.GetAttributeName(i))
2344 if(
"MDMoveEnergy"==tag.GetAttributeName(i))
2348 if(
"LogLikelihoodScale"==tag.GetAttributeName(i))
2356 if((
"Molecule"==tagg.GetName())&&tagg.IsEndTag())
2360 VFN_DEBUG_EXIT(
"Molecule::XMLInput():"<<this->
GetName(),5)
2363 if(
"Quaternion"==tagg.GetName())
2365 mQuat.XMLInput(is,tagg);
2367 if(
"Atom"==tagg.GetName())
2370 mvpAtom.back()->XMLInput(is,tagg);
2372 if(
"Bond"==tagg.GetName())
2374 this->
AddBond(this->GetAtom(0),this->GetAtom(1),1.5,.01,.05,1.,
false);
2375 mvpBond.back()->XMLInput(is,tagg);
2377 if(
"BondAngle"==tagg.GetName())
2379 this->
AddBondAngle(this->GetAtom(0),this->GetAtom(1),this->GetAtom(2),1.5,.01,.05,
false);
2382 if(
"DihedralAngle"==tagg.GetName())
2385 this->GetAtom(2),this->GetAtom(3),1.5,.01,.05,
false);
2388 if(
"RigidGroup"==tagg.GetName())
2391 for(
unsigned int i=0;i<tagg.GetNbAttribute();i++)
2392 if(tagg.GetAttributeName(i).rfind(
"Atom", 0)==0)
2393 s.insert(&(this->GetAtom(tagg.GetAttributeValue(i))));
2396 if(
"CenterAtom"==tagg.GetName())
2399 for(
unsigned int i=0;i<tagg.GetNbAttribute();i++)
2400 if(
"Name"==tagg.GetAttributeName(i))
2401 this->
SetCenterAtom(this->GetAtom(tagg.GetAttributeValue(i)));
2403 if(
"Option"==tagg.GetName())
2405 for(
unsigned int i=0;i<tagg.GetNbAttribute();i++)
2406 if(
"Name"==tagg.GetAttributeName(i))
2410 if(
"Par"==tagg.GetName())
2412 for(
unsigned int i=0;i<tagg.GetNbAttribute();i++)
2414 if(
"Name"==tagg.GetAttributeName(i))
2416 if(
"x"==tagg.GetAttributeValue(i))
2421 if(
"y"==tagg.GetAttributeValue(i))
2426 if(
"z"==tagg.GetAttributeValue(i))
2431 if(
"Occup"==tagg.GetAttributeValue(i))
2440 VFN_DEBUG_EXIT(
"Molecule::XMLInput()",5)
2457 TAU_PROFILE(
"Molecule::BeginOptimization()",
"void (bool,bool)",TAU_DEFAULT);
2459 if( (!mIsSelfOptimizing)
2464 (*fpObjCrystInformUser)(
"Optimizing initial conformation of Molecule:"+this->
GetName());
2466 (*fpObjCrystInformUser)(
"Finished optimizing initial conformation of Molecule:"+this->
GetName());
2473 clockConf=mClockAtomList;
2474 if(clockConf<mClockBondList) clockConf=mClockBondList;
2475 if(clockConf<mClockBondAngleList) clockConf=mClockBondAngleList;
2476 if(clockConf<mClockDihedralAngleList) clockConf=mClockDihedralAngleList;
2477 if(clockConf<mClockRigidGroup) clockConf=mClockRigidGroup;
2478 if(clockConf<mClockAtomScattPow) clockConf=mClockAtomScattPow;
2480 clockMode=mClockConnectivityTable;
2481 if(clockMode<mClockRingList) clockMode=mClockRingList;
2482 if(clockMode<mClockRotorGroup) clockMode=mClockRotorGroup;
2483 if(clockMode<mClockFlipGroup) clockMode=mClockFlipGroup;
2484 if(clockMode<mClockStretchModeBondLength) clockMode=mClockStretchModeBondLength;
2485 if(clockMode<mClockStretchModeBondAngle) clockMode=mClockStretchModeBondAngle;
2486 if(clockMode<mClockStretchModeTorsion) clockMode=mClockStretchModeTorsion;
2487 if(clockMode<mClockStretchModeTwist) clockMode=mClockStretchModeTwist;
2488 if(clockMode<mClockMDAtomGroup) clockMode=mClockMDAtomGroup;
2491 if( (!mIsSelfOptimizing) && (clockMode<clockConf))
2525 for(vector<MolAtom*>::iterator at=this->GetAtomList().begin();at!=this->GetAtomList().end();++at)
2527 this->
GetPar(&((*at)->X())).SetIsFixed(
true);
2528 this->
GetPar(&((*at)->Y())).SetIsFixed(
true);
2529 this->
GetPar(&((*at)->Z())).SetIsFixed(
true);
2534 for(vector<MolAtom*>::iterator at=this->GetAtomList().begin();at!=this->GetAtomList().end();++at)
2536 this->
GetPar(&((*at)->X())).SetIsFixed(
false);
2537 this->
GetPar(&((*at)->Y())).SetIsFixed(
false);
2538 this->
GetPar(&((*at)->Z())).SetIsFixed(
false);
2541 #ifdef RIGID_BODY_STRICT_EXPERIMENTAL
2550 (*pos)->mQuat.Q0()=1;
2551 (*pos)->mQuat.Q1()=0;
2552 (*pos)->mQuat.Q2()=0;
2553 (*pos)->mQuat.Q3()=0;
2554 (*pos)->mvIdx.clear();
2555 for(set<MolAtom *>::iterator at=(*pos)->begin();at!=(*pos)->end();++at)
2557 this->
GetPar(&((*at)->X())).SetIsFixed(
true);
2558 this->
GetPar(&((*at)->Y())).SetIsFixed(
true);
2559 this->
GetPar(&((*at)->Z())).SetIsFixed(
true);
2561 if(&(this->GetAtom(i))==*at)
2563 (*pos)->mvIdx.insert(i);
2571 mRandomConformChangeNbTest=0;
2572 mRandomConformChangeNbAccept=0;
2573 mRandomConformChangeTemp=1.;
2585 #ifdef RIGID_BODY_STRICT_EXPERIMENTAL
2589 for(set<MolAtom *>::iterator at=(*pos)->begin();at!=(*pos)->end();++at)
2591 this->
GetPar(&((*at)->X())).SetIsFixed(
false);
2592 this->
GetPar(&((*at)->Y())).SetIsFixed(
false);
2593 this->
GetPar(&((*at)->Z())).SetIsFixed(
false);
2605 TAU_PROFILE(
"Molecule::RandomizeConfiguration()",
"void ()",TAU_DEFAULT);
2606 VFN_DEBUG_ENTRY(
"Molecule::RandomizeConfiguration()",4)
2608 if( (!mIsSelfOptimizing)
2612 (*fpObjCrystInformUser)(
"Optimizing initial conformation of Molecule:"+this->
GetName());
2614 (*fpObjCrystInformUser)(
"Finished optimizing initial conformation of Molecule:"+this->
GetName());
2643 const REAL angle=(REAL)rand()*2.*M_PI/(REAL)RAND_MAX;
2645 pos->mvRotatedAtomList,angle);
2651 const REAL angle=(REAL)rand()*2.*M_PI/(REAL)RAND_MAX;
2653 pos->mvRotatedAtomList,angle);
2659 const REAL angle=(REAL)rand()*2.*M_PI/(REAL)RAND_MAX;
2661 pos->mvRotatedAtomList,angle);
2664 for(list<StretchModeTorsion>::const_iterator
2668 const REAL amp=2*M_PI*rand()/(REAL)RAND_MAX;
2675 map<MolAtom*,XYZ> v0;
2676 for(vector<MolAtom*>::iterator at=this->GetAtomList().begin();at!=this->GetAtomList().end();++at)
2677 v0[*at]=
XYZ(rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5);
2680 +this->GetBondAngleList().size()
2681 +this->GetDihedralAngleList().size());
2682 map<RigidGroup*,std::pair<XYZ,XYZ> > vr;
2684 this->GetBondList(),
2685 this->GetBondAngleList(),
2686 this->GetDihedralAngleList(),
2693 #ifdef RIGID_BODY_STRICT_EXPERIMENTAL
2701 (*pos)->mQuat.Q0()=1;
2702 (*pos)->mQuat.Q1()=0;
2703 (*pos)->mQuat.Q2()=0;
2704 (*pos)->mQuat.Q3()=0;
2709 const REAL amp=M_PI/RAND_MAX;
2711 ((2.*(REAL)rand()-(REAL)RAND_MAX)*amp,
2712 (REAL)rand(),(REAL)rand(),(REAL)rand());
2714 mClockOrientation.
Click();
2716 VFN_DEBUG_EXIT(
"Molecule::RandomizeConfiguration()",4)
2723 if(mIsSelfOptimizing)
2727 VFN_DEBUG_EXIT(
"Molecule::GlobalOptRandomMove()",4)
2730 TAU_PROFILE(
"Molecule::GlobalOptRandomMove()",
"void (REAL,RefParType*)",TAU_DEFAULT);
2731 TAU_PROFILE_TIMER(timer1,
"Molecule::GlobalOptRandomMove 1",
"", TAU_FIELD);
2732 TAU_PROFILE_TIMER(timer2,
"Molecule::GlobalOptRandomMove 2",
"", TAU_FIELD);
2733 TAU_PROFILE_TIMER(timer3,
"Molecule::GlobalOptRandomMove 3",
"", TAU_FIELD);
2734 TAU_PROFILE_TIMER(timer4,
"Molecule::GlobalOptRandomMove 4",
"", TAU_FIELD);
2735 TAU_PROFILE_TIMER(timer5,
"Molecule::GlobalOptRandomMove 5",
"", TAU_FIELD);
2736 VFN_DEBUG_ENTRY(
"Molecule::GlobalOptRandomMove()",4)
2743 &&(mFlipModel.GetChoice()==0)
2746 &&(((rand()%100)==0)))
2751 const unsigned long i=rand() %
mvFlipGroup.size();
2752 list<FlipGroup>::iterator pos=
mvFlipGroup.begin();
2753 for(
unsigned long j=0;j<i;++j)++pos;
2756 static unsigned long ctflip1=0,ctflip2=0;
2757 if(pos->mvRotatedChainList.begin()->first==pos->mpAtom0)
2759 cout <<
"TRYING: Flip group from atom "
2760 <<pos->mpAtom0->GetName()<<
",exchanging bonds with "
2761 <<pos->mpAtom1->GetName()<<
" and "
2762 <<pos->mpAtom2->GetName()<<
", resulting in a 180deg rotation of atoms : ";
2763 for(set<MolAtom*>::iterator pos1=pos->mvRotatedChainList.begin()->second.begin();
2764 pos1!=pos->mvRotatedChainList.begin()->second.end();++pos1)
2765 cout<<(*pos1)->GetName()<<
" ";
2769 cout <<
"TRYING: Flip group with respect to: "
2770 <<pos->mpAtom1->GetName()<<
"-"
2771 <<pos->mpAtom0->GetName()<<
"-"
2772 <<pos->mpAtom2->GetName()<<
" : ";
2773 for(list<pair<
const MolAtom *,set<MolAtom*> > >::const_iterator
2774 chain=pos->mvRotatedChainList.begin();
2775 chain!=pos->mvRotatedChainList.end();++chain)
2777 cout<<
" -"<<chain->first->GetName()<<
":";
2778 for(set<MolAtom*>::const_iterator pos1=chain->second.begin();
2779 pos1!=chain->second.end();++pos1)
2780 cout<<(*pos1)->GetName()<<
" ";
2792 cout<<
" FLIP ACCEPTED ("<<float(ctflip2)/ctflip1*100<<
"% accepted)"<<endl;
2802 TAU_PROFILE_START(timer1);
2805 static const REAL amp=mBaseRotationAmplitude/(REAL)RAND_MAX;
2809 ((2.*(REAL)rand()-(REAL)RAND_MAX)*amp*mutationAmplitude*mult,
2810 (REAL)rand(),(REAL)rand(),(REAL)rand());
2812 mClockOrientation.
Click();
2822 TAU_PROFILE_STOP(timer1);
2841 if((*at)->GetX()<xmin) xmin=(*at)->GetX();
2842 if((*at)->GetX()>xmax) xmax=(*at)->GetX();
2843 if((*at)->GetY()<ymin) ymin=(*at)->GetY();
2844 if((*at)->GetY()>ymax) ymax=(*at)->GetY();
2845 if((*at)->GetZ()<zmin) zmin=(*at)->GetZ();
2846 if((*at)->GetZ()>zmax) zmax=(*at)->GetZ();
2849 REAL dx=(xmax-xmin)/5.,dy=(ymax-ymin)/5.,dz=(zmax-zmin)/5.;
2853 const REAL xc=xmin+rand()/(REAL)RAND_MAX*(xmax-xmin);
2854 const REAL yc=ymin+rand()/(REAL)RAND_MAX*(ymax-ymin);
2855 const REAL zc=zmin+rand()/(REAL)RAND_MAX*(zmax-zmin);
2856 map<MolAtom*,XYZ> v0;
2857 const REAL ax=-4.*log(2.)/(dx*dx);
2858 const REAL ay=-4.*log(2.)/(dy*dy);
2859 const REAL az=-4.*log(2.)/(dz*dz);
2861 v0[*at]=
XYZ(exp(ax*((*at)->GetX()-xc)*((*at)->GetX()-xc)),
2862 exp(ay*((*at)->GetY()-yc)*((*at)->GetY()-yc)),
2863 exp(az*((*at)->GetZ()-zc)*((*at)->GetZ()-zc)));
2866 map<MolAtom*,XYZ> v0;
2869 std::map<MolAtom*,unsigned long> pushedAtoms;
2870 unsigned long idx=rand()%v0.size();
2872 for(
unsigned int i=0;i<idx;i++) at0++;
2873 const REAL xc=(*at0)->GetX();
2874 const REAL yc=(*at0)->GetY();
2875 const REAL zc=(*at0)->GetZ();
2881 ux=REAL(rand()-RAND_MAX/2);
2882 uy=REAL(rand()-RAND_MAX/2);
2883 uz=REAL(rand()-RAND_MAX/2);
2884 n=sqrt(ux*ux+uy*uy+uz*uz);
2886 ux=ux/n;uy=uy/n;uz=uz/n;
2887 const REAL a=-4.*log(2.)/(2*2);
2889 for(map<MolAtom*,unsigned long>::iterator at=pushedAtoms.begin() ;at!=pushedAtoms.end();++at)
2890 v0[at->first]=
XYZ(ux*exp(a*(at->first->GetX()-xc)*(at->first->GetX()-xc)),
2891 uy*exp(a*(at->first->GetY()-yc)*(at->first->GetY()-yc)),
2892 uz*exp(a*(at->first->GetZ()-zc)*(at->first->GetZ()-zc)));
2894 for(map<MolAtom*, unsigned long>::iterator at=pushedAtoms.begin() ;at!=pushedAtoms.end();++at)
2895 v0[at->first]=
XYZ((at->first->GetX()-xc)*ux*exp(a*(at->first->GetX()-xc)*(at->first->GetX()-xc)),
2896 (at->first->GetY()-yc)*uy*exp(a*(at->first->GetY()-yc)*(at->first->GetY()-yc)),
2897 (at->first->GetZ()-zc)*uz*exp(a*(at->first->GetZ()-zc)*(at->first->GetZ()-zc)));
2902 +this->GetBondAngleList().size()
2903 +this->GetDihedralAngleList().size());
2904 map<RigidGroup*,std::pair<XYZ,XYZ> > vr;
2906 this->GetBondList(),
2907 this->GetBondAngleList(),
2908 this->GetDihedralAngleList(),
2916 for(
unsigned int i=0;i<n;++i)++pos;
2917 map<MolAtom*,XYZ> v0;
2918 for(set<MolAtom*>::iterator at=pos->mvpAtom.begin();at!=pos->mvpAtom.end();++at)
2919 v0[*at]=
XYZ(rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5);
2922 +pos->mvpBondAngle.size()
2923 +pos->mvpDihedralAngle.size());
2924 map<RigidGroup*,std::pair<XYZ,XYZ> > vr;
2925 float nrjMult=1.0+mutationAmplitude*0.2;
2926 if((rand()%20)==0) nrjMult=4.0;
2930 pos->mvpDihedralAngle,
2938 for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
2950 (*mode)->CalcDeriv();
2952 for(map<const MolBond*,REAL>::const_iterator pos=(*mode)->mvpBrokenBond.begin();
2953 pos!=(*mode)->mvpBrokenBond.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
2954 for(map<const MolBondAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenBondAngle.begin();
2955 pos!=(*mode)->mvpBrokenBondAngle.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
2956 for(map<const MolDihedralAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenDihedralAngle.begin();
2957 pos!=(*mode)->mvpBrokenDihedralAngle.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
2959 REAL change=(2.*(REAL)rand()-(REAL)RAND_MAX)/(REAL)RAND_MAX;
2962 if((*mode)->mLLKDeriv>0)
2964 change -= 0.3*sqrt(llk);
2965 if(change<-1) change=-1;
2969 change += 0.3*sqrt(llk);
2970 if(change>1) change=1;
2972 (*mode)->Print(cout);
2973 change *= mutationAmplitude * (*mode)->mBaseAmplitude;
2974 cout <<
" Change="<<change<<
" (dLLK= "<<(*mode)->mLLKDeriv<<
"), llk= "<<llk<<
" ?->"<<llk+(*mode)->mLLKDeriv*change<<endl;
2976 (*mode)->Stretch(change);
2979 for(map<const MolBond*,REAL>::const_iterator pos=(*mode)->mvpBrokenBond.begin();
2980 pos!=(*mode)->mvpBrokenBond.end();++pos)
2982 cout<<
" "<<pos->first->GetName()<<
", llk= "<<pos->first->GetLogLikelihood(
false,
false)
2983 <<
" ?->"<<pos->first->GetLogLikelihood(
false,
false)+pos->first->GetDeriv((*mode)->mDerivXYZ,
true)*change
2984 <<
"? (deriv="<<pos->first->GetDeriv((*mode)->mDerivXYZ)<<
", "<<pos->first->GetDeriv((*mode)->mDerivXYZ,
true);
2985 cout<<
") ->" <<pos->first->GetLogLikelihood()<<endl;
2986 llk+=pos->first->GetLogLikelihood(
false,
false);
2988 for(map<const MolBondAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenBondAngle.begin();
2989 pos!=(*mode)->mvpBrokenBondAngle.end();++pos)
2991 cout<<
" "<<pos->first->GetName()<<
", llk= "<<pos->first->GetLogLikelihood(
false,
false)
2992 <<
" ?->"<<pos->first->GetLogLikelihood(
false,
false)+pos->first->GetDeriv((*mode)->mDerivXYZ,
true)*change
2993 <<
"? (deriv="<<pos->first->GetDeriv((*mode)->mDerivXYZ)<<
", "<<pos->first->GetDeriv((*mode)->mDerivXYZ,
true);
2994 cout<<
") ->" <<pos->first->GetLogLikelihood()<<endl;
2995 llk+=pos->first->GetLogLikelihood(
false,
false);
2997 for(map<const MolDihedralAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenDihedralAngle.begin();
2998 pos!=(*mode)->mvpBrokenDihedralAngle.end();++pos)
3000 cout<<
" "<<pos->first->GetName()<<
", llk= "<<pos->first->GetLogLikelihood(
false,
false)
3001 <<
" ?->"<<pos->first->GetLogLikelihood(
false,
false)+pos->first->GetDeriv((*mode)->mDerivXYZ,
true)*change
3002 <<
"? (deriv="<<pos->first->GetDeriv((*mode)->mDerivXYZ)<<
", "<<pos->first->GetDeriv((*mode)->mDerivXYZ,
true);
3003 cout<<
") ->" <<pos->first->GetLogLikelihood()<<endl;
3004 llk+=pos->first->GetLogLikelihood(
false,
false);
3006 cout <<
" -> "<<llk<<endl;
3011 TAU_PROFILE_START(timer2);
3015 if((rand()%2)==0) (*mode)->RandomStretch(mutationAmplitude);
3017 TAU_PROFILE_STOP(timer2);
3023 TAU_PROFILE_START(timer3);
3024 for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
3030 TAU_PROFILE_STOP(timer3);
3032 TAU_PROFILE_START(timer4);
3040 (*mode)->CalcDeriv();
3042 for(map<const MolBond*,REAL>::const_iterator pos=(*mode)->mvpBrokenBond.begin();
3043 pos!=(*mode)->mvpBrokenBond.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
3044 for(map<const MolBondAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenBondAngle.begin();
3045 pos!=(*mode)->mvpBrokenBondAngle.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
3046 for(map<const MolDihedralAngle*,REAL>::const_iterator pos=(*mode)->mvpBrokenDihedralAngle.begin();
3047 pos!=(*mode)->mvpBrokenDihedralAngle.end();++pos) llk+=pos->first->GetLogLikelihood(
false,
false);
3048 REAL change=(2.*(REAL)rand()-(REAL)RAND_MAX)/(REAL)RAND_MAX;
3050 if((*mode)->mLLKDeriv>0)
3053 if(change<-1) change=-1;
3058 if(change>1) change=1;
3060 if(mutationAmplitude<0.5) change *= mutationAmplitude * (*mode)->mBaseAmplitude;
3061 else change *= 0.5 * (*mode)->mBaseAmplitude;
3062 (*mode)->Stretch(change);
3070 TAU_PROFILE_STOP(timer4);
3079 map<MolAtom*,XYZ> v0;
3080 for(set<MolAtom*>::iterator at=pos->mvpAtom.begin();at!=pos->mvpAtom.end();++at)
3081 v0[*at]=
XYZ(rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5,rand()/(REAL)RAND_MAX+0.5);
3083 const REAL nrj0=20*(pos->mvpBond.size()+pos->mvpBondAngle.size()+pos->mvpDihedralAngle.size());
3084 map<RigidGroup*,std::pair<XYZ,XYZ> > vr;
3086 (
const vector<MolBond*>) (pos->mvpBond),
3087 (
const vector<MolBondAngle*>) (pos->mvpBondAngle),
3088 (
const vector<MolDihedralAngle*>) (pos->mvpDihedralAngle),
3097 mClockLogLikelihood.
Click();
3104 REAL x0=0,y0=0,z0=0;
3105 for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
3114 for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
3122 VFN_DEBUG_EXIT(
"Molecule::GlobalOptRandomMove()",4)
3127 if( (mClockLogLikelihood>mClockAtomList)
3128 &&(mClockLogLikelihood>mClockBondList)
3129 &&(mClockLogLikelihood>mClockBondAngleList)
3130 &&(mClockLogLikelihood>mClockDihedralAngleList)
3131 &&(mClockLogLikelihood>mClockAtomPosition)
3133 TAU_PROFILE(
"Molecule::GetLogLikelihood()",
"REAL ()",TAU_DEFAULT);
3135 mClockLogLikelihood.
Click();
3148 for(vector<MolBond*>::const_iterator pos=this->GetBondList().begin();pos!=this->GetBondList().end();++pos)
3149 *p++=(*pos)->GetLength();
3150 for(vector<MolBondAngle*>::const_iterator pos=this->GetBondAngleList().begin();pos!=this->GetBondAngleList().end();++pos)
3151 *p++=(*pos)->GetAngle();
3152 for(vector<MolDihedralAngle*>::const_iterator pos=this->GetDihedralAngleList().begin();pos!=this->GetDihedralAngleList().end();++pos)
3153 *p++=(*pos)->GetAngle();
3161 for(vector<MolBond*>::const_iterator pos=this->GetBondList().begin();pos!=this->GetBondList().end();++pos)
3162 *p++=(*pos)->GetLength0();
3163 for(vector<MolBondAngle*>::const_iterator pos=this->GetBondAngleList().begin();pos!=this->GetBondAngleList().end();++pos)
3164 *p++=(*pos)->GetAngle0();
3165 for(vector<MolDihedralAngle*>::const_iterator pos=this->GetDihedralAngleList().begin();pos!=this->GetDihedralAngleList().end();++pos)
3166 *p++=(*pos)->GetAngle0();
3175 for(vector<MolBond*>::const_iterator pos=this->GetBondList().begin();pos!=this->GetBondList().end();++pos)
3176 *p++=1/((*pos)->GetLengthSigma()* (*pos)->GetLengthSigma()+1e-6);
3177 for(vector<MolBondAngle*>::const_iterator pos=this->GetBondAngleList().begin();pos!=this->GetBondAngleList().end();++pos)
3178 *p++=1/((*pos)->GetAngleSigma()* (*pos)->GetAngleSigma()+1e-6);
3179 for(vector<MolDihedralAngle*>::const_iterator pos=this->GetDihedralAngleList().begin();pos!=this->GetDihedralAngleList().end();++pos)
3180 *p++=1/((*pos)->GetAngleSigma()* (*pos)->GetAngleSigma()+1e-6);
3195 cout<<
"Molecule::TagNewBestConfig()"<<endl;
3197 vector<MolBond*>::const_iterator pos;
3200 cout<<
"BondLength="<<(*pos)->GetLength();
3201 (*pos)->XMLOutput(cout);
3205 vector<MolBondAngle*>::const_iterator pos;
3208 cout<<
"BondAngle="<<(*pos)->GetAngle();
3209 (*pos)->XMLOutput(cout);
3213 vector<MolDihedralAngle*>::const_iterator pos;
3216 cout<<
"DihedralAngle="<<(*pos)->GetAngle();
3217 (*pos)->XMLOutput(cout);
3224 for(set<MolAtom*>::iterator at1=pos->mvpAtom.begin();at1!=pos->mvpAtom.end();++at1)
3226 sprintf(buf,
"%5s : ",(*at1)->GetName().c_str());
3228 for(set<MolAtom*>::iterator at2=at1;at2!=pos->mvpAtom.end();++at2)
3230 if(at1==at2)
continue;
3231 sprintf(buf,
"%5s(%6.3f),",(*at2)->GetName().c_str(),
GetBondLength(**at1,**at2));
3244 VFN_DEBUG_ENTRY(
"Molecule::GetScatteringComponentList()",3)
3246 VFN_DEBUG_EXIT(
"Molecule::GetScatteringComponentList()",3)
3258 VFN_DEBUG_ENTRY(
"Molecule::POVRayDescription()",3)
3259 const REAL xMin=options.
mXmin;
const REAL xMax=options.mXmax;
3260 const REAL yMin=options.mYmin;
const REAL yMax=options.mYmax;
3261 const REAL zMin=options.mZmin;
const REAL zMax=options.mZmax;
3264 VFN_DEBUG_EXIT(
"Molecule::POVRayDescription():No atom to display !",4)
3273 os <<
"// Description of Molecule :" << this->
GetName() <<endl;
3274 vector<CrystMatrix_REAL> vXYZCoords;
3283 vXYZCoords.push_back(this->
GetCrystal().GetSpaceGroup().
3284 GetAllSymmetrics(x0,y0,z0,
false,
false,
false));
3287 CrystMatrix_int translate(27,3);
3288 translate= -1,-1,-1,
3317 CrystVector_REAL xSave,ySave,zSave;
3318 const int nbSymmetrics=vXYZCoords[0].rows();
3320 for(
int i=0;i<nbSymmetrics;i++)
3322 VFN_DEBUG_ENTRY(
"Molecule::POVRayDescription():Symmetric#"<<i,3)
3323 for(
unsigned int j=0;j<
mvpAtom.size();j++)
3325 x(j)=vXYZCoords[j](i,0);
3326 y(j)=vXYZCoords[j](i,1);
3327 z(j)=vXYZCoords[j](i,2);
3333 x(0) = fmod((
float) x(0),(
float)1);
if(x(0)<0) x(0)+=1.;
3334 y(0) = fmod((
float) y(0),(
float)1);
if(y(0)<0) y(0)+=1.;
3335 z(0) = fmod((
float) z(0),(
float)1);
if(z(0)<0) z(0)+=1.;
3339 for(
unsigned int j=1;j<
mvpAtom.size();j++)
3349 for(
int j=0;j<translate.rows();j++)
3351 x += translate(j,0);
3352 y += translate(j,1);
3353 z += translate(j,2);
3358 if( ((x.min()<xMax) && (x.max()>xMin))
3359 &&((y.min()<yMax) && (y.max()>yMin))
3360 &&((z.min()<zMax) && (z.max()>zMin)))
3362 os<<
" //Symetric#"<<++ct<<endl;
3363 for(
unsigned int k=0;k<
mvpAtom.size();k++)
3365 isinside(k)=((x(k)>=xMin) && (x(k)<=xMax)) && ((y(k)>=yMin) && (y(k)<=yMax)) && ((z(k)>=zMin) && (z(k)<=zMax));
3366 if(isinside(k)) borderdist(k)=0;
3370 if(xMin>x(k)) borderdist(k)+=(xMin-x(k))*aa*(xMin-x(k))*aa;
3371 if(yMin>y(k)) borderdist(k)+=(yMin-y(k))*bb*(yMin-y(k))*bb;
3372 if(zMin>z(k)) borderdist(k)+=(zMin-z(k))*cc*(zMin-z(k))*cc;
3373 if(xMax<x(k)) borderdist(k)+=(xMax-x(k))*aa*(xMax-x(k))*aa;
3374 if(yMax<y(k)) borderdist(k)+=(yMax-y(k))*bb*(yMax-y(k))*bb;
3375 if(zMax<z(k)) borderdist(k)+=(zMax-z(k))*cc*(zMax-z(k))*cc;
3376 borderdist(k)=sqrt(borderdist(k));
3382 if((
mvpAtom[k]->IsDummy()) || (fout<0.001))
continue;
3383 if(options.
mShowHydrogens==
false && (
mvpAtom[k]->GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5))
continue;
3408 os <<
" ObjCrystAtom("
3412 <<
mvpAtom[k]->GetScatteringPower().GetRadius()/3.0<<
","
3413 <<
"colour_"+
mvpAtom[k]->GetScatteringPower().GetName()<<
","
3417 for(
unsigned int k=0;k<
mvpBond.size();k++)
3419 if( (
mvpBond[k]->GetAtom1().IsDummy())
3420 ||(
mvpBond[k]->GetAtom2().IsDummy()) )
continue;
3421 if(options.
mShowHydrogens==
false && ( (
mvpBond[k]->GetAtom1().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)
3422 ||(
mvpBond[k]->GetAtom2().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)))
continue;
3423 unsigned long n1,n2;
3425 for(n1=0;n1<
mvpAtom.size();n1++)
3427 for(n2=0;n2<
mvpAtom.size();n2++)
3430 if((isinside(n1)==
false) || (isinside(n2)==
false))
3432 if(fout<0.001)
continue;
3433 REAL x1=x(n1),y1=y(n1),z1=z(n1),
3434 x2=x(n2),y2=y(n2),z2=z(n2);
3435 REAL dx=x2-x1,dy=y2-y1,dz=z2-z1;
3436 const REAL r=sqrt(abs(dx*dx+dy*dy+dz*dz))+1e-6;
3437 const REAL r1=
mvpAtom[n1]->GetScatteringPower().GetRadius()/3.0;
3438 const REAL r2=
mvpAtom[n2]->GetScatteringPower().GetRadius()/3.0;
3439 x1+=dx/r*r1*sqrt(abs(1-0.1/r1));
3440 y1+=dy/r*r1*sqrt(abs(1-0.1/r1));
3441 z1+=dz/r*r1*sqrt(abs(1-0.1/r1));
3442 x2-=dx/r*r2*sqrt(abs(1-0.1/r2));
3443 y2-=dy/r*r2*sqrt(abs(1-0.1/r2));
3444 z2-=dz/r*r2*sqrt(abs(1-0.1/r2));
3446 os <<
" ObjCrystBond("
3447 <<x1<<
","<<y1<<
","<<z1<<
","
3448 <<x2<<
","<<y2<<
","<<z2<<
","
3450 if(
mvpBond[k]->IsFreeTorsion()) os<<
"colour_freebond,";
3451 else os<<
"colour_nonfreebond,";
3452 os<<f<<
","<<fout<<
")"<<endl;
3459 VFN_DEBUG_EXIT(
"Molecule::POVRayDescription():Symmetric#"<<i,3)
3461 VFN_DEBUG_EXIT(
"Molecule::POVRayDescription()",3)
3466 void Molecule::GLInitDisplayList(
const bool onlyIndependentAtoms,
3467 const REAL xMin,
const REAL xMax,
3468 const REAL yMin,
const REAL yMax,
3469 const REAL zMin,
const REAL zMax,
3470 const bool displayEnantiomer,
3471 const bool displayNames,
3472 const bool hideHydrogens,
3473 const REAL fadeDistance,
3474 const bool fullMoleculeInLimits)
const
3476 VFN_DEBUG_ENTRY(
"Molecule::GLInitDisplayList()",3)
3479 VFN_DEBUG_EXIT(
"Molecule::GLInitDisplayList():No atom to display !",4)
3484 if(
mvpAtom.size()>200) large=true;
3486 if(displayEnantiomer==true) en=-1;
3493 const REAL aa=this->
GetCrystal().GetLatticePar(0);
3494 const REAL bb=this->
GetCrystal().GetLatticePar(1);
3495 const REAL cc=this->
GetCrystal().GetLatticePar(2);
3497 const GLfloat colour0[] = {0.0f, 0.0f, 0.0f, 0.0f};
3498 glMaterialfv(GL_FRONT, GL_SPECULAR, colour0);
3499 glMaterialfv(GL_FRONT, GL_EMISSION, colour0);
3500 glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
3501 glPolygonMode(GL_FRONT, GL_FILL);
3503 GLUquadricObj* pQuadric = gluNewQuadric();
3505 if(
true==onlyIndependentAtoms)
3509 vector<MolAtom*>::const_iterator pos;
3513 if((*pos)->IsDummy())
continue;
3514 if(hideHydrogens && ((*pos)->GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5))
continue;
3515 const float r=(*pos)->GetScatteringPower().GetColourRGB()[0];
3516 const float g=(*pos)->GetScatteringPower().GetColourRGB()[1];
3517 const float b=(*pos)->GetScatteringPower().GetColourRGB()[2];
3518 const float f=(*pos)->GetOccupancy()*this->
GetOccupancy();
3522 GLfloat colourChar [] = {1.0, 1.0, 1.0, 1.0};
3523 GLfloat colourCharRing [] = {1.0, 1.0, 0.8, 1.0};
3524 if((r>0.8)&&(g>0.8)&&(b>0.8))
3526 colourChar[0] = 0.5;
3527 colourChar[1] = 0.5;
3528 colourChar[2] = 0.5;
3530 if((*pos)->IsInRing())
3531 glMaterialfv(GL_FRONT, GL_EMISSION, colourCharRing);
3533 glMaterialfv(GL_FRONT, GL_EMISSION, colourChar);
3534 glMaterialfv(GL_FRONT, GL_SHININESS, colour0);
3535 glTranslatef((*pos)->X()*en+xc, (*pos)->Y()+yc, (*pos)->Z()+zc);
3536 crystGLPrint((*pos)->GetName());
3540 const GLfloat colourAtom [] = {r, g, b, f};
3541 glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, colourAtom);
3542 glTranslatef((*pos)->X()*en+xc, (*pos)->Y()+yc, (*pos)->Z()+zc);
3543 gluSphere(pQuadric,(*pos)->GetScatteringPower().GetRadius()/3.,20,20);
3550 VFN_DEBUG_ENTRY(
"Molecule::GLInitDisplayList():Show all symmetrics",3)
3552 map<const MolAtom*,
unsigned long> rix;
3555 for(vector<MolAtom*>::const_iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
3558 vector<CrystMatrix_REAL> vXYZCoords;
3567 vXYZCoords.push_back(this->
GetCrystal().GetSpaceGroup().
3568 GetAllSymmetrics(x0,y0,z0,
false,
false,
false));
3571 CrystMatrix_int translate(27,3);
3572 translate= -1,-1,-1,
3601 CrystVector_REAL xSave,ySave,zSave;
3602 const int nbSymmetrics=vXYZCoords[0].rows();
3603 for(
int i=0;i<nbSymmetrics;i++)
3605 VFN_DEBUG_ENTRY(
"Molecule::GLInitDisplayList():Symmetric#"<<i,3)
3606 for(
unsigned int j=0;j<
mvpAtom.size();j++)
3608 x(j)=vXYZCoords[j](i,0);
3609 y(j)=vXYZCoords[j](i,1);
3610 z(j)=vXYZCoords[j](i,2);
3616 x(0) = fmod((
float) x(0),(
float)1);
if(x(0)<0) x(0)+=1.;
3617 y(0) = fmod((
float) y(0),(
float)1);
if(y(0)<0) y(0)+=1.;
3618 z(0) = fmod((
float) z(0),(
float)1);
if(z(0)<0) z(0)+=1.;
3622 for(
unsigned int j=1;j<
mvpAtom.size();j++)
3632 for(
int j=0;j<translate.rows();j++)
3634 x += translate(j,0);
3635 y += translate(j,1);
3636 z += translate(j,2);
3639 const bool molcenter_isinside =( ((x.sum()/x.size())>=xMin) && ((x.sum()/x.size())<=xMax)
3640 && ((y.sum()/y.size())>=yMin) && ((y.sum()/y.size())<=yMax)
3641 && ((z.sum()/z.size())>=zMin) && ((z.sum()/z.size())<=zMax));
3642 if( ((x.min()<(xMax+fadeDistance/aa)) && (x.max()>(xMin-fadeDistance/aa)))
3643 &&((y.min()<(yMax+fadeDistance/bb)) && (y.max()>(yMin-fadeDistance/bb)))
3644 &&((z.min()<(zMax+fadeDistance/cc)) && (z.max()>(zMin-fadeDistance/cc))))
3646 for(
unsigned int k=0;k<
mvpAtom.size();k++)
3648 isinside(k)=((x(k)>=xMin) && (x(k)<=xMax)) && ((y(k)>=yMin) && (y(k)<=yMax)) && ((z(k)>=zMin) && (z(k)<=zMax));
3649 if(isinside(k)) borderdist(k)=0;
3653 if(xMin>x(k)) borderdist(k)+=(xMin-x(k))*aa*(xMin-x(k))*aa;
3654 if(yMin>y(k)) borderdist(k)+=(yMin-y(k))*bb*(yMin-y(k))*bb;
3655 if(zMin>z(k)) borderdist(k)+=(zMin-z(k))*cc*(zMin-z(k))*cc;
3656 if(xMax<x(k)) borderdist(k)+=(xMax-x(k))*aa*(xMax-x(k))*aa;
3657 if(yMax<y(k)) borderdist(k)+=(yMax-y(k))*bb*(yMax-y(k))*bb;
3658 if(zMax<z(k)) borderdist(k)+=(zMax-z(k))*cc*(zMax-z(k))*cc;
3659 borderdist(k)=sqrt(borderdist(k));
3662 if((isinside(k)==
false) && ((molcenter_isinside==
false) || (fullMoleculeInLimits==
false)))
3664 if ((fadeDistance == 0) || borderdist(k)>fadeDistance) fout = 0;
3669 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);
3670 VFN_DEBUG_MESSAGE(ch,4)
3673 if(
mvpAtom[k]->IsDummy())
continue;
3674 if(hideHydrogens && (
mvpAtom[k]->GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5))
continue;
3675 if(fout<0.01)
continue;
3677 const float r=
mvpAtom[k]->GetScatteringPower().GetColourRGB()[0];
3678 const float g=
mvpAtom[k]->GetScatteringPower().GetColourRGB()[1];
3679 const float b=
mvpAtom[k]->GetScatteringPower().GetColourRGB()[2];
3683 if((fout>0.99) || molcenter_isinside)
3685 GLfloat colourChar [] = {1.0, 1.0, 1.0, f*fout};
3686 GLfloat colourCharRing [] = {1.0, 1.0, 0.8, f*fout};
3687 if((r>0.8)&&(g>0.8)&&(b>0.8))
3689 colourChar[0] = 0.5;
3690 colourChar[1] = 0.5;
3691 colourChar[2] = 0.5;
3694 glMaterialfv(GL_FRONT, GL_EMISSION, colourCharRing);
3696 glMaterialfv(GL_FRONT, GL_EMISSION, colourChar);
3697 glRasterPos3f(x(k)*en, y(k), z(k));
3705 const GLfloat colourAtom [] = {r*fout, g*fout, b*fout, f*fout};
3706 glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,colourAtom);
3707 glTranslatef(x(k)*en, y(k), z(k));
3709 mvpAtom[k]->GetScatteringPower().GetRadius()/3.,20,20);
3713 const GLfloat colourAtom [] = {r*fout, g*fout, b*fout, f*fout};
3714 glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,colourAtom);
3715 glTranslatef(x(k)*en, y(k), z(k));
3716 gluSphere(pQuadric,.15,10,10);
3721 if(displayNames==
false)
3725 for(
unsigned int k=0;k<
mvpBond.size();k++)
3727 if( (
mvpBond[k]->GetAtom1().IsDummy())
3728 ||(
mvpBond[k]->GetAtom2().IsDummy()) )
continue;
3729 if(hideHydrogens && ( (
mvpBond[k]->GetAtom1().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)
3730 ||(
mvpBond[k]->GetAtom2().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)))
continue;
3731 const unsigned long n1=rix[&(
mvpBond[k]->GetAtom1())],
3732 n2=rix[&(
mvpBond[k]->GetAtom2())];
3734 if(((isinside(n1)==
false) || (isinside(n2)==
false)) && ((molcenter_isinside==
false) || (fullMoleculeInLimits==
false)))
3736 if((fadeDistance==0) || ((borderdist(n1)+borderdist(n2))/2)>fadeDistance) fout = 0;
3737 else fout*=(fadeDistance-(borderdist(n1)+borderdist(n2))/2)/fadeDistance*
3740 if(fout<0.01)
continue;
3742 const float r1=
mvpBond[k]->GetAtom1().GetScatteringPower().GetColourRGB()[0];
3743 const float g1=
mvpBond[k]->GetAtom1().GetScatteringPower().GetColourRGB()[1];
3744 const float b1=
mvpBond[k]->GetAtom1().GetScatteringPower().GetColourRGB()[2];
3746 const float r2=
mvpBond[k]->GetAtom2().GetScatteringPower().GetColourRGB()[0];
3747 const float g2=
mvpBond[k]->GetAtom2().GetScatteringPower().GetColourRGB()[1];
3748 const float b2=
mvpBond[k]->GetAtom2().GetScatteringPower().GetColourRGB()[2];
3750 const GLfloat colourAtom1 [] = {r1, g1, b1, f1*fout};
3751 const GLfloat colourAtom2 [] = {r2, g2, b2, f2*fout};
3754 glBegin(GL_LINE_STRIP);
3755 glMaterialfv(GL_FRONT, GL_SPECULAR, colourAtom1);
3756 glMaterialfv(GL_FRONT, GL_EMISSION, colourAtom1);
3757 glMaterialfv(GL_FRONT, GL_SHININESS, colourAtom1);
3758 glVertex3f(x(n1)*en,y(n1),z(n1));
3759 glVertex3f((x(n1)+x(n2))/2*en,(y(n1)+y(n2))/2,(z(n1)+z(n2))/2);
3760 glMaterialfv(GL_FRONT, GL_SPECULAR, colourAtom2);
3761 glMaterialfv(GL_FRONT, GL_EMISSION, colourAtom2);
3762 glMaterialfv(GL_FRONT, GL_SHININESS, colourAtom2);
3763 glVertex3f(x(n2)*en,y(n2),z(n2));
3767 const REAL height= sqrt(abs( (x(n2)-x(n1))*(x(n2)-x(n1))
3768 +(y(n2)-y(n1))*(y(n2)-y(n1))
3769 +(z(n2)-z(n1))*(z(n2)-z(n1))));
3771 glTranslatef(x(n1)*en, y(n1), z(n1));
3772 GLUquadricObj *quadobj = gluNewQuadric();
3773 glRotatef(180,(x(n2)-x(n1))*en,y(n2)-y(n1),z(n2)-z(n1)+height);
3775 glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colourAtom1);
3776 gluCylinder(quadobj,.1,.1,height/2,10,1 );
3777 gluDeleteQuadric(quadobj);
3779 glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colourAtom2);
3780 GLUquadricObj *quadobj2 = gluNewQuadric();
3781 glTranslatef(0, 0, height/2);
3782 gluCylinder(quadobj2,.1,.1,height/2,10,1 );
3783 gluDeleteQuadric(quadobj2);
3790 for(
unsigned int k=0;k<
mvpBond.size();k++)
3792 if( (
mvpBond[k]->GetAtom1().IsDummy())
3793 ||(
mvpBond[k]->GetAtom2().IsDummy()) )
continue;
3794 if(hideHydrogens && ( (
mvpBond[k]->GetAtom1().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)
3795 ||(
mvpBond[k]->GetAtom2().GetScatteringPower().GetForwardScatteringFactor(RAD_XRAY)<1.5)))
continue;
3796 const unsigned long n1=rix[&(
mvpBond[k]->GetAtom1())],
3797 n2=rix[&(
mvpBond[k]->GetAtom2())];
3800 if(((isinside(n1)==
false) || (isinside(n2)==
false)) && ((molcenter_isinside==
false) || (fullMoleculeInLimits==
false)))
3802 if ((fadeDistance == 0) || ((borderdist(n1) + borderdist(n2)) / 2)>fadeDistance)
3805 fout=(fadeDistance-(borderdist(n1)+borderdist(n2))/2)/fadeDistance*
3809 if(fout<0.01)
continue;
3811 bool isRigidGroup=
false;
3814 if( ((*pos)->find(&(
mvpBond[k]->GetAtom1()))!=(*pos)->end())
3815 &&((*pos)->find(&(
mvpBond[k]->GetAtom2()))!=(*pos)->end()) )
3822 const GLfloat colour_bondnonfree[]= { 0.2f*fout, .2f*fout, .2f*fout, f*fout };
3823 const GLfloat colour_bondrigid[]= { 0.5f*fout, .3f*fout, .3f*fout, f*fout };
3824 const GLfloat colour_bondfree[]= { 0.8f*fout, .8f*fout, .8f*fout, f*fout };
3826 glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colour_bondrigid);
3829 if(
mvpBond[k]->IsFreeTorsion())
3830 glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colour_bondfree);
3832 glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,colour_bondnonfree);
3835 REAL x1=x(n1),y1=y(n1),z1=z(n1),
3836 x2=x(n2),y2=y(n2),z2=z(n2);
3837 REAL dx=x2-x1,dy=y2-y1,dz=z2-z1;
3838 const REAL r=sqrt(abs(dx*dx+dy*dy+dz*dz))+1e-6;
3839 const REAL r1=
mvpAtom[n1]->GetScatteringPower().GetRadius()/3.0;
3840 const REAL r2=
mvpAtom[n2]->GetScatteringPower().GetRadius()/3.0;
3841 x1+=dx/r*r1*sqrt(abs(1-0.1/r1));
3842 y1+=dy/r*r1*sqrt(abs(1-0.1/r1));
3843 z1+=dz/r*r1*sqrt(abs(1-0.1/r1));
3844 x2-=dx/r*r2*sqrt(abs(1-0.1/r2));
3845 y2-=dy/r*r2*sqrt(abs(1-0.1/r2));
3846 z2-=dz/r*r2*sqrt(abs(1-0.1/r2));
3849 glTranslatef(x1*en, y1, z1);
3850 GLUquadricObj *quadobj = gluNewQuadric();
3852 const REAL height= sqrt(abs( (x2-x1)*(x2-x1)
3855 glRotatef(180,(x2-x1)*en,y2-y1,z2-z1+height);
3856 gluCylinder(quadobj,.1,.1,height,10,1 );
3857 gluDeleteQuadric(quadobj);
3867 VFN_DEBUG_EXIT(
"Molecule::GLInitDisplayList():Symmetric#"<<i,3)
3869 VFN_DEBUG_EXIT(
"Molecule::GLInitDisplayList():Show all symmetrics",3)
3871 gluDeleteQuadric(pQuadric);
3872 VFN_DEBUG_EXIT(
"Molecule::GLInitDisplayList()",3)
3878 const bool updateDisplay)
3880 VFN_DEBUG_ENTRY(
"Molecule::AddAtom():"<<name,5)
3883 string thename=name;
3884 if(thename==
string(
""))
3887 if(pPow!=0) sprintf(buf,
"%s_%s_%lu",this->
GetName().c_str(),pPow->
GetName().c_str(),
mvpAtom.size()+1);
3888 else sprintf(buf,
"%s_X_%lu",this->
GetName().c_str(),
mvpAtom.size()+1);
3892 mClockAtomPosition.
Click();
3893 mClockAtomScattPow.
Click();
3897 gpRefParTypeScattConformX,
3898 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
3906 gpRefParTypeScattConformY,
3907 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
3915 gpRefParTypeScattConformZ,
3916 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
3925 if(molnameasformula || (this->
GetName().size()==0))
3929 VFN_DEBUG_EXIT(
"Molecule::AddAtom()",5)
3934 VFN_DEBUG_ENTRY(
"Molecule::RemoveAtom():"<<atom.GetName(),6)
3935 vector<MolAtom*>::iterator pos=find(
mvpAtom.begin(),
mvpAtom.end(),&atom);
3939 +
" is not in this Molecule:"+this->GetName());
3947 for(vector<MolBond*>::iterator posb=
mvpBond.begin();posb!=
mvpBond.end();)
3949 if( (&atom==&((*posb)->GetAtom1())) || (&atom==&((*posb)->GetAtom2())) )
3957 if( (&atom==&((*posb)->GetAtom1())) || (&atom==&((*posb)->GetAtom2()))
3958 ||(&atom==&((*posb)->GetAtom3())))
3967 if( (&atom==&((*posb)->GetAtom1())) || (&atom==&((*posb)->GetAtom2()))
3968 ||(&atom==&((*posb)->GetAtom3())) || (&atom==&((*posb)->GetAtom4())))
3972 mClockAtomList.
Click();
3976 if(del)
delete *pos;
3980 if(molnameasformula || (this->
GetName().size()==0))
3984 VFN_DEBUG_EXIT(
"Molecule::RemoveAtom()",6)
3989 const REAL length,
const REAL sigma,
const REAL delta,
3990 const REAL bondOrder,
3991 const bool updateDisplay)
3993 VFN_DEBUG_ENTRY(
"Molecule::AddBond()",5)
3994 mvpBond.push_back(
new MolBond(atom1,atom2,length,sigma,delta,*
this,bondOrder));
3996 mClockBondList.
Click();
3998 VFN_DEBUG_EXIT(
"Molecule::AddBond()",5)
4003 VFN_DEBUG_ENTRY(
"Molecule::RemoveBond():"<<bond.
GetName(),6)
4004 vector<MolBond*>::iterator pos=find(
mvpBond.begin(),
mvpBond.end(),&bond);
4008 +
"-"+bond.GetAtom2().GetName()
4009 +
" is not in this Molecule:"+this->GetName());
4012 mClockBondList.
Click();
4013 if(del)
delete *pos;
4016 VFN_DEBUG_EXIT(
"Molecule::RemoveBond():",6)
4022 for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
4024 if( ((&((*pos)->GetAtom1())==&at1)&&(&((*pos)->GetAtom2())==&at2))
4025 ||((&((*pos)->GetAtom1())==&at2)&&(&((*pos)->GetAtom2())==&at1)))
4032 for(vector<MolBond*>::iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
4034 if( ((&((*pos)->GetAtom1())==&at1)&&(&((*pos)->GetAtom2())==&at2))
4035 ||((&((*pos)->GetAtom1())==&at2)&&(&((*pos)->GetAtom2())==&at1)))
4042 const REAL angle,
const REAL sigma,
const REAL delta,
4043 const bool updateDisplay)
4045 VFN_DEBUG_ENTRY(
"Molecule::AddBondAngle()",5)
4048 mClockBondAngleList.
Click();
4050 VFN_DEBUG_EXIT(
"Molecule::AddBondAngle()",5)
4055 VFN_DEBUG_ENTRY(
"Molecule::RemoveBondAngle():"<<angle.GetName(),6)
4059 throw ObjCrystException(
"Molecule::RemoveBondAngle():"+angle.GetAtom1().GetName()
4060 +
"-"+angle.GetAtom2().GetName()+
"-"+angle.GetAtom3().GetName()
4061 +
" is not in this Molecule:"+this->GetName());
4064 mClockBondAngleList.
Click();
4065 if(del)
delete *pos;
4068 VFN_DEBUG_EXIT(
"Molecule::RemoveBondAngle():",6)
4076 for(vector<MolBondAngle*>::const_iterator pos=
mvpBondAngle.begin();
4079 if( (&((*pos)->GetAtom2())==&at2)
4080 &&( ((&((*pos)->GetAtom1())==&at1)&&(&((*pos)->GetAtom3())==&at3))
4081 ||((&((*pos)->GetAtom1())==&at3)&&(&((*pos)->GetAtom3())==&at1))))
4089 const REAL angle,
const REAL sigma,
const REAL delta,
4090 const bool updateDisplay)
4092 VFN_DEBUG_ENTRY(
"Molecule::AddDihedralAngle()",5)
4094 angle,sigma,delta,*
this));
4096 mClockDihedralAngleList.
Click();
4098 VFN_DEBUG_EXIT(
"Molecule::AddDihedralAngle()",5)
4103 VFN_DEBUG_ENTRY(
"Molecule::RemoveDihedralAngle():"<<angle.GetName(),6)
4108 throw ObjCrystException(
"Molecule::RemoveDihedralAngle():"+angle.GetAtom1().GetName()
4109 +
"-"+angle.GetAtom2().GetName()+
"-"+angle.GetAtom3().GetName()
4110 +
"-"+angle.GetAtom4().GetName()
4111 +
" is not in this Molecule:"+this->GetName());
4114 mClockDihedralAngleList.
Click();
4115 if(del)
delete *pos;
4118 VFN_DEBUG_ENTRY(
"Molecule::RemoveDihedralAngle():",6)
4126 for(vector<MolDihedralAngle*>::const_iterator pos=
mvpDihedralAngle.begin();
4129 if( ( ((&((*pos)->GetAtom1())==&at1)&&(&((*pos)->GetAtom2())==&at2))
4130 &&((&((*pos)->GetAtom3())==&at3)&&(&((*pos)->GetAtom4())==&at4)))
4131 ||( ((&((*pos)->GetAtom4())==&at1)&&(&((*pos)->GetAtom3())==&at2))
4132 &&((&((*pos)->GetAtom2())==&at3)&&(&((*pos)->GetAtom1())==&at4))))
4139 const bool updateDisplay)
4142 #ifdef RIGID_BODY_STRICT_EXPERIMENTAL
4154 sprintf(buf,
"RigidGroup%d_x",i);
4156 gpRefParTypeScattConformX,
4157 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
4164 sprintf(buf,
"RigidGroup%d_y",i);
4166 gpRefParTypeScattConformY,
4167 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
4174 sprintf(buf,
"RigidGroup%d_z",i);
4176 gpRefParTypeScattConformZ,
4177 REFPAR_DERIV_STEP_ABSOLUTE,
false,
false,
true,
false,1.,1.);
4184 sprintf(buf,
"RigidGroup%d_Q1",i);
4186 gpRefParTypeScattConform,
4187 REFPAR_DERIV_STEP_ABSOLUTE,
true,
false,
true,
false,1.,1.);
4194 sprintf(buf,
"RigidGroup%d_Q2",i);
4196 gpRefParTypeScattConform,
4197 REFPAR_DERIV_STEP_ABSOLUTE,
true,
false,
true,
false,1.,1.);
4204 sprintf(buf,
"RigidGroup%d_Q3",i);
4206 gpRefParTypeScattConform,
4207 REFPAR_DERIV_STEP_ABSOLUTE,
true,
false,
true,
false,1.,1.);
4214 mClockRigidGroup.
Click();
4222 #ifdef RIGID_BODY_STRICT_EXPERIMENTAL
4235 if(del)
delete *pos;
4243 const MolAtom &Molecule::GetAtom(
unsigned int i)
const{
return *
mvpAtom[i];}
4245 MolAtom &Molecule::GetAtom(
const string &name){
return **(this->
FindAtom(name));}
4247 const MolAtom &Molecule::GetAtom(
const string &name)
const{
return **(this->
FindAtom(name));}
4251 VFN_DEBUG_ENTRY(
"Molecule::OptimizeConformation()",5)
4255 ANNEALING_EXPONENTIAL,10,.1);
4258 mIsSelfOptimizing=
true;
4259 globalOptObj.
Optimize(nb,
false,stopCost);
4260 mIsSelfOptimizing=
false;
4262 mClockFlipGroup.
Reset();
4263 mClockRotorGroup.
Reset();
4264 VFN_DEBUG_EXIT(
"Molecule::OptimizeConformation()",5)
4270 for(
unsigned i = 0; i < nbStep; ++i)
4274 map<MolAtom*,XYZ> grad;
4275 for(vector<MolAtom*>::iterator pos=this->GetAtomList().begin();
4276 pos!=this->GetAtomList().end();++pos) grad[*pos]=
XYZ(0,0,0);
4279 for(vector<MolBond*>::iterator pos=this->GetBondList().begin();
4280 pos!=this->GetBondList().end();++pos) (*pos)->CalcGradient(grad);
4282 for(vector<MolBondAngle*>::iterator pos=this->GetBondAngleList().begin();
4283 pos!=this->GetBondAngleList().end();++pos) (*pos)->CalcGradient(grad);
4285 for(vector<MolDihedralAngle*>::iterator pos=this->GetDihedralAngleList().begin();
4286 pos!=this->GetDihedralAngleList().end();++pos) (*pos)->CalcGradient(grad);
4290 for(map<MolAtom*,XYZ>::const_iterator pos=grad.begin();pos!=grad.end();++pos)
4293 sprintf(buf,
"%10s Grad LLK= (%8.3f %8.3f %8.3f)",
4294 pos->first->GetName().c_str(),pos->second.x,pos->second.y,pos->second.z);
4300 for(map<MolAtom*,XYZ>::const_iterator pos=grad.begin();pos!=grad.end();++pos)
4302 if(abs(pos->second.x)>f) f=abs(pos->second.x);
4303 if(abs(pos->second.y)>f) f=abs(pos->second.y);
4304 if(abs(pos->second.z)>f) f=abs(pos->second.z);
4306 if(f>1e-6) f=maxStep/f;
4313 if((*pos)->size()==0)
continue;
4314 REAL dx=0,dy=0,dz=0;
4315 for(set<MolAtom *>::const_iterator at=(*pos)->begin();at!=(*pos)->end();++at)
4324 for(set<MolAtom *>::const_iterator at=(*pos)->begin();at!=(*pos)->end();++at)
4332 for(map<MolAtom*,XYZ>::const_iterator pos=grad.begin();pos!=grad.end();++pos)
4334 pos->first->SetX(pos->first->GetX()-pos->second.x*f);
4335 pos->first->SetY(pos->first->GetY()-pos->second.y*f);
4336 pos->first->SetZ(pos->first->GetZ()-pos->second.z*f);
4344 const vector<MolBond*> &vb,
const vector<MolBondAngle*> &va,
4345 const vector<MolDihedralAngle*> &vd,
4346 map<
RigidGroup*,std::pair<XYZ,XYZ> > &vr, REAL nrj0)
4348 const vector<MolBond*> *pvb=&vb;
4349 const vector<MolBondAngle*> *pva=&va;
4350 const vector<MolDihedralAngle*> *pvd=&vd;
4351 map<RigidGroup*,std::pair<XYZ,XYZ> > *pvr=&vr;
4352 if((pvb->size()==0)&&(pva->size()==0)&&(pvd->size()==0))
4354 pvb=&(this->GetBondList());
4355 pva=&(this->GetBondAngleList());
4356 pvd=&(this->GetDihedralAngleList());
4358 (*pvr)[*pos]=make_pair(
XYZ(0,0,0),
XYZ(0,0,0));
4364 REAL e_v,e_k,v_r=1.0;
4365 for(
unsigned i = 0; i < nbStep; ++i)
4368 map<MolAtom*,XYZ> grad;
4369 for(map<MolAtom*,XYZ>::iterator pos=v0.begin();pos!=v0.end();++pos) grad[pos->first]=
XYZ(0,0,0);
4373 for(vector<MolBond*>::const_iterator pos=pvb->begin();pos!=pvb->end();++pos)
4375 (*pos)->CalcGradient(grad);
4376 e_v+=(*pos)->GetLogLikelihood(
false,
false);
4379 for(vector<MolBondAngle*>::const_iterator pos=pva->begin();pos!=pva->end();++pos)
4381 (*pos)->CalcGradient(grad);
4382 e_v+=(*pos)->GetLogLikelihood(
false,
false);
4385 for(vector<MolDihedralAngle*>::const_iterator pos=pvd->begin();pos!=pvd->end();++pos)
4387 (*pos)->CalcGradient(grad);
4388 e_v+=(*pos)->GetLogLikelihood(
false,
false);
4393 for(map<MolAtom*,XYZ>::const_iterator pos=v0.begin();pos!=v0.end();++pos)
4394 e_k += 0.5*m*(pos->second.x*pos->second.x + pos->second.y*pos->second.y + pos->second.z*pos->second.z);
4396 if(nrj0==0) nrj0=e_k+e_v;
4400 const REAL de=e_k+e_v-nrj0;
4401 if(de<e_k) v_r=sqrt((e_k-de)/e_k);
4407 sprintf(buf,
"(i) LLK + Ek = %10.3f + %10.3f =%10.3f (nrj0=%10.3f)",e_v,e_k,e_v+e_k,nrj0);
4411 for(map<MolAtom*,XYZ>::iterator pos=v0.begin();pos!=v0.end();++pos)
4413 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)",
4414 pos->first->GetName().c_str(),
4415 pos->first->GetX(),pos->first->GetY(),pos->first->GetZ(),
4416 v0[*pos].x ,v0[*pos].y ,v0[*pos].z,
4417 -grad[*pos].x*im,-grad[*pos].y*im,-grad[*pos].z*im);
4423 for(map<MolAtom*,XYZ>::const_iterator pos=v0.begin();pos!=v0.end();++pos)
4425 const XYZ *pa=&(grad[pos->first]);
4426 pos->first->SetX(pos->first->GetX()+pos->second.x*dt*v_r-0.5*im*dt*dt*pa->x);
4427 pos->first->SetY(pos->first->GetY()+pos->second.y*dt*v_r-0.5*im*dt*dt*pa->y);
4428 pos->first->SetZ(pos->first->GetZ()+pos->second.z*dt*v_r-0.5*im*dt*dt*pa->z);
4431 for(map<MolAtom*,XYZ>::iterator pos=v0.begin();pos!=v0.end();++pos)
4433 const XYZ *pa=&(grad[pos->first]);
4434 pos->second.x = v_r*pos->second.x - pa->x*dt*im;
4435 pos->second.y = v_r*pos->second.y - pa->y*dt*im;
4436 pos->second.z = v_r*pos->second.z - pa->z*dt*im;
4441 const vector<MolAtom*>& Molecule::GetAtomList()
const{
return mvpAtom;}
4442 const vector<MolBond*>& Molecule::GetBondList()
const{
return mvpBond;}
4443 const vector<MolBondAngle*>& Molecule::GetBondAngleList()
const{
return mvpBondAngle;}
4444 const vector<MolDihedralAngle*>& Molecule::GetDihedralAngleList()
const{
return mvpDihedralAngle;}
4446 vector<MolAtom*>& Molecule::GetAtomList(){
return mvpAtom;}
4447 vector<MolBond*>& Molecule::GetBondList(){
return mvpBond;}
4448 vector<MolBondAngle*>& Molecule::GetBondAngleList(){
return mvpBondAngle;}
4449 vector<MolDihedralAngle*>& Molecule::GetDihedralAngleList(){
return mvpDihedralAngle;}
4455 const list<StretchModeBondLength>& Molecule::GetStretchModeBondLengthList()
const{
return mvStretchModeBondLength;}
4456 const list<StretchModeBondAngle>& Molecule::GetStretchModeBondAngleList()
const{
return mvStretchModeBondAngle;}
4457 const list<StretchModeTorsion>& Molecule::GetStretchModeTorsionList()
const{
return mvStretchModeTorsion;}
4463 const set<MolAtom *> &atoms,
const REAL angle,
4464 const bool keepCenter)
4466 const REAL vx=at2.X()-at1.X();
4467 const REAL vy=at2.Y()-at1.Y();
4468 const REAL vz=at2.Z()-at1.Z();
4472 const set<MolAtom *> &atoms,
const REAL angle,
4473 const bool keepCenter)
4475 TAU_PROFILE(
"Molecule::RotateAtomGroup(MolAtom&,vx,vy,vz,...)",
"void (...)",TAU_DEFAULT);
4476 if(atoms.size()==0)
return;
4477 const REAL x0=at.X();
4478 const REAL y0=at.Y();
4479 const REAL z0=at.Z();
4481 if((fabs(vx)+fabs(vy)+fabs(vz))<1e-6)
return;
4482 REAL dx=0.,dy=0.,dz=0.;
4483 bool keepc=keepCenter;
4486 ||(this->GetPar(
mXYZ.data()+1).IsFixed())
4487 ||(this->GetPar(
mXYZ.data()+2).IsFixed())) keepc=
false;
4489 const REAL ca=cos(angle),sa=sin(angle);
4490 const REAL ca1=1-ca;
4491 const REAL vnorm=1/sqrt(vx*vx+vy*vy+vz*vz);
4492 const REAL ux=vx*vnorm,uy=vy*vnorm,uz=vz*vnorm;
4493 const REAL m00=ca+ux*ux*ca1;
4494 const REAL m01=ux*uy*ca1-uz*sa;
4495 const REAL m02=ux*uz*ca1+uy*sa;
4496 const REAL m10=uy*ux*ca1+uz*sa;
4497 const REAL m11=ca+uy*uy*ca1;
4498 const REAL m12=uy*uz*ca1-ux*sa;
4499 const REAL m20=uz*ux*ca1-uy*sa;
4500 const REAL m21=uz*uy*ca1+ux*sa;
4501 const REAL m22=ca+uz*uz*ca1;
4502 for(set<MolAtom *>::const_iterator pos=atoms.begin();pos!=atoms.end();++pos)
4510 const REAL x=(*pos)->X() - x0;
4511 const REAL y=(*pos)->Y() - y0;
4512 const REAL z=(*pos)->Z() - z0;
4514 (*pos)->X() = m00*x+m01*y+m02*z+x0;
4515 (*pos)->Y() = m10*x+m11*y+m12*z+y0;
4516 (*pos)->Z() = m20*x+m21*y+m22*z+z0;
4526 for(set<MolAtom *>::const_iterator pos=atoms.begin();pos!=atoms.end();++pos)
4537 quat.
RotateVector((*pos)->X(),(*pos)->Y(),(*pos)->Z());
4561 mClockAtomPosition.
Click();
4565 const REAL dx,
const REAL dy,
const REAL dz,
4566 const bool keepCenter)
4568 for(set<MolAtom *>::const_iterator pos=atoms.begin();pos!=atoms.end();++pos)
4574 bool keepc=keepCenter;
4577 ||(this->GetPar(
mXYZ.data()+1).IsFixed())
4578 ||(this->GetPar(
mXYZ.data()+2).IsFixed())) keepc=
false;
4581 const REAL r= (REAL)(atoms.size())/(REAL)(this->
GetNbComponent());
4582 REAL dxc=dx*r,dyc=dy*r,dzc=dz*r;
4588 mClockAtomPosition.
Click();
4593 VFN_DEBUG_ENTRY(
"Molecule::RestraintExport()",5)
4594 os<<
"BondName, IdealLength, Length, log(likelihood)"<<endl;
4595 for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
4596 os <<(*pos)->GetName()
4597 <<
", "<<(*pos)->GetLength0()
4598 <<
", "<<(*pos)->GetLength()
4599 <<
", "<<(*pos)->GetLogLikelihood()<<endl;
4600 os<<
"BondAngle, IdealAngle, Angle, log(likelihood)"<<endl;
4601 for(vector<MolBondAngle*>::const_iterator pos=
mvpBondAngle.begin();
4603 os <<(*pos)->GetName()
4604 <<
", "<<(*pos)->Angle0()*180/M_PI
4605 <<
", "<<(*pos)->GetAngle()*180/M_PI
4606 <<
", "<<(*pos)->GetLogLikelihood()<<endl;
4607 os<<
"DihedralAngle, IdealAngle, Angle, log(likelihood)"<<endl;
4608 for(vector<MolDihedralAngle*>::const_iterator pos=
mvpDihedralAngle.begin();
4610 os <<(*pos)->GetName()
4611 <<
", "<<(*pos)->Angle0()*180/M_PI
4612 <<
", "<<(*pos)->GetAngle()*180/M_PI
4613 <<
", "<<(*pos)->GetLogLikelihood()<<endl;
4614 VFN_DEBUG_EXIT(
"Molecule::RestraintExport()",5)
4618 VFN_DEBUG_ENTRY(
"Molecule::RestraintStatus()",5)
4619 for(vector<MolBond*>::const_iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
4620 cout <<
"Bond "<<(*pos)->GetName()
4621 <<
", IdealLength="<<(*pos)->GetLength0()
4622 <<
", Length="<<(*pos)->GetLength()
4623 <<
", log(likelihood)="<<(*pos)->GetLogLikelihood()<<endl;
4624 for(vector<MolBondAngle*>::const_iterator pos=
mvpBondAngle.begin();
4626 cout <<
"Bond Angle "<<(*pos)->GetName()
4627 <<
", IdealAngle="<<(*pos)->Angle0()*180/M_PI
4628 <<
", Angle="<<(*pos)->GetAngle()*180/M_PI
4629 <<
", log(likelihood)="<<(*pos)->GetLogLikelihood()<<endl;
4630 for(vector<MolDihedralAngle*>::const_iterator pos=
mvpDihedralAngle.begin();
4632 cout <<
"Dihedral Angle "<<(*pos)->GetName()
4633 <<
", IdealAngle="<<(*pos)->Angle0()*180/M_PI
4634 <<
", Angle="<<(*pos)->GetAngle()*180/M_PI
4635 <<
", log(likelihood)="<<(*pos)->GetLogLikelihood()<<endl;
4636 VFN_DEBUG_EXIT(
"Molecule::RestraintStatus()",5)
4660 for(vector<MolBond*>::iterator bond=
mvpBond.begin();bond!=
mvpBond.end();++bond)
4662 MolAtom* at2=&((*bond)->GetAtom1());
4663 MolAtom* at3=&((*bond)->GetAtom2());
4668 if(*c2==at3)
continue;
4669 if(
GetBondAngle(**c2,*at2,*at3)<(10 *DEG2RAD))
continue;
4670 if(
GetBondAngle(**c2,*at2,*at3)>(180*DEG2RAD))
continue;
4675 if((*c3==at2)||(*c3==*c2))
continue;
4676 if(
GetBondAngle(*at2,*at3,**c3)<(10 *DEG2RAD))
continue;
4677 if(
GetBondAngle(*at2,*at3,**c3)>(180*DEG2RAD))
continue;
4692 static const REAL svGaussianIntX[51]
4693 ={0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2,
4694 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2, 2.3,
4695 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3. , 3.1, 3.2, 3.3, 3.4,
4696 3.5, 3.6, 3.7, 3.8, 3.9, 4. , 4.1, 4.2, 4.3, 4.4, 4.5,
4697 4.6, 4.7, 4.8, 4.9, 5. };
4701 static const REAL svGaussianIntY[51]
4702 ={0.00000000e+00, 1.00285768e-01, 2.02498389e-01, 3.08782899e-01,
4703 4.21537858e-01, 5.43577677e-01, 6.78339751e-01, 8.30161516e-01,
4704 1.00466359e+00, 1.20929269e+00, 1.45410564e+00, 1.75292008e+00,
4705 2.12502806e+00, 2.59778353e+00, 3.21056320e+00, 4.02091257e+00,
4706 5.11421527e+00, 6.61911896e+00, 8.73249677e+00, 1.17604260e+01,
4707 1.61864569e+01, 2.27870547e+01, 3.28297946e+01, 4.84189023e+01,
4708 7.31071567e+01, 1.12996748e+02, 1.78751757e+02, 2.89337246e+02,
4709 4.79080996e+02, 8.11232960e+02, 1.40443983e+03, 2.48531490e+03,
4710 4.49461494e+03, 8.30539633e+03, 1.56790580e+04, 3.02354014e+04,
4711 5.95525189e+04, 1.19793234e+05, 2.46080284e+05, 5.16182008e+05,
4712 1.10556243e+06, 2.41765307e+06, 5.39775929e+06, 1.23033271e+07,
4713 2.86288375e+07, 6.80050408e+07, 1.64899866e+08, 4.08157783e+08,
4714 1.03122228e+09, 2.65938808e+09, 7.00012860e+09};
4718 static const CubicSpline sInvGaussianInt(svGaussianIntY,svGaussianIntX,51);
4720 REAL FlatGaussianProba(
const REAL x,
const REAL sigma,
const REAL delta)
4722 if(abs(x)<delta)
return 1.0;
4723 const REAL z=(abs(x)-delta)/sigma;
4727 REAL FlatGaussianIntegral(
const REAL x1,
const REAL x2,
const REAL sigma,
const REAL delta)
4729 static const REAL SPI2=0.88622692545275794;
4730 if((x1<=-delta)&&(x2<=-delta))
return SPI2*sigma*(erf((x2+delta)/sigma)-erf((x1+delta)/sigma));
4731 if((x1<=-delta)&&(x2<= delta))
return (x2+delta)-SPI2*sigma*erf((x1+delta)/sigma);
4732 if(x1<=-delta)
return 2*delta+SPI2*sigma*(erf((x2-delta)/sigma)-erf((x1+delta)/sigma));
4733 if((x1<= delta)&&(x2<= delta))
return x2-x1;
4734 if(x1<= delta)
return SPI2*sigma*erf((x2-delta)/sigma)+(delta-x1);
4735 return SPI2*sigma*(erf((x2-delta)/sigma)-erf((x1-delta)/sigma));
4738 REAL FlatLorentzianProba(
const REAL x,
const REAL sigma,
const REAL delta)
4740 if(abs(x)<delta)
return 1.0;
4741 const REAL z=(abs(x)-delta)/sigma;
4745 REAL FlatLorentzianIntegral(
const REAL x1,
const REAL x2,
const REAL sigma,
const REAL delta)
4747 if((x1<=-delta)&&(x2<=-delta))
return atan((x2+delta)/sigma)-atan((x1+delta)/sigma);
4748 if((x1<=-delta)&&(x2<= delta))
return (x2+delta)-atan((x1+delta)/sigma);
4749 if(x1<=-delta)
return 2*delta+atan((x2-delta)/sigma)-atan((x1+delta)/sigma);
4750 if((x1<= delta)&&(x2<= delta))
return x2-x1;
4751 if(x1<= delta)
return atan((x2-delta)/sigma)+(delta-x1);
4752 return atan((x2-delta)/sigma)-atan((x1-delta)/sigma);
4764 REAL r=(REAL)rand()/(REAL)RAND_MAX;
4767 REAL x=x0+amplitude*(2*r-1.0);
4768 if(x> delta)x= delta;
4769 if(x<-delta)x=-delta;
4781 Pmin=FlatLorentzianProba((x0+xmin)/2,sigma,delta);
4782 Pmax=FlatLorentzianProba((x0+xmax)/2,sigma,delta);
4787 for(
int i=0;i<5;i++)
4789 const REAL r=Pmax/Pmin;
4790 xmax=x0+r*amplitude;
4791 Pmax=FlatLorentzianProba((x0+xmax)/2,sigma,delta);
4797 for(
int i=0;i<5;i++)
4799 const REAL r=Pmin/Pmax;
4800 xmin=x0-r*amplitude;
4801 Pmin=FlatLorentzianProba((x0+xmin)/2,sigma,delta);
4814 REAL d=(abs(x0)-delta)/sigma;
4816 d=(1-amplitude*d/2)/(1+amplitude*d/2);
4819 xmin=x0-amplitude*2/(1+d);
4820 xmax=x0+amplitude*2*d/(1+d);
4824 xmax=x0+amplitude*2/(1+d);
4825 xmin=x0-amplitude*2*d/(1+d);
4834 REAL ymin=(abs(xmin)-delta)/sigma;
4836 REAL ymax=(abs(xmax)-delta)/sigma;
4838 const REAL y=ymin+(ymax-ymin)*r;
4839 return -tan(y)*sigma-delta;
4846 REAL p0=atan((abs(xmin)-delta)/sigma)*sigma;
4851 REAL ymin=(abs(xmin)-delta)/sigma;
4853 const REAL y=ymin*(REAL)rand()/(REAL)RAND_MAX;
4854 return -delta-tan(y)*sigma;
4858 return -delta+(REAL)rand()/(REAL)RAND_MAX*(xmax+delta);
4863 REAL p0=atan((abs(xmin)-delta)/sigma)*sigma;
4865 REAL p2=atan((xmax-delta)/sigma)*sigma;
4866 const REAL n=p0+p1+p2;
4869 REAL ymin=(abs(xmin)-delta)/sigma;
4871 const REAL y=ymin*(REAL)rand()/(REAL)RAND_MAX;
4872 const REAL x=-delta-tan(y)*sigma;
4877 const REAL x=-delta+(REAL)rand()/(REAL)RAND_MAX*2*delta;
4881 REAL ymax=(xmax-delta)/sigma;
4883 const REAL y=ymax*(REAL)rand()/(REAL)RAND_MAX;
4884 const REAL x=delta+tan(y)*sigma;
4892 return xmin+r*(xmax-xmin);
4895 const REAL p0=delta-xmin;
4896 const REAL p1=atan((xmax-delta)/sigma)*sigma;
4899 return xmin+(REAL)rand()/(REAL)RAND_MAX*(delta-xmin);
4902 REAL ymax=(xmax-delta)/sigma;
4904 const REAL y=ymax*(REAL)rand()/(REAL)RAND_MAX;
4905 return delta+tan(y)*sigma;
4908 REAL ymin=(xmin-delta)/sigma;
4910 REAL ymax=(xmax-delta)/sigma;
4912 const REAL y=ymin+(ymax-ymin)*r;
4913 return tan(y)*sigma+delta;
4916 void TestLorentzianBiasedRandomMove()
4919 REAL x=0,sigma=0.1,delta=0.5,amplitude=0.05;
4922 for(
long i=0;i<400000;i++)
4955 const bool respectRestraint)
4960 const REAL l=sqrt(dx*dx+dy*dy+dz*dz+1e-7);
4962 if(l<1e-6)
return change;
4963 if(respectRestraint && mode.
mpBond!=0)
4965 const REAL d0=l-mode.
mpBond->GetLength0();
4966 const REAL sigma=mode.
mpBond->GetLengthSigma();
4967 const REAL delta=mode.
mpBond->GetLengthDelta();
4968 const REAL max=delta+sigma*5.0;
4971 REAL d1=d0+(REAL)(2*rand()-RAND_MAX)/(REAL)RAND_MAX*amplitude*0.1;
4972 if(d1> delta)d1= delta;
4973 if(d1<-delta)d1=-delta;
4977 if((d0+change)>max) change=max-d0;
4978 else if((d0+change)<(-max)) change=-max-d0;
4982 cout<<
"BOND LENGTH change("<<change<<
"):"
4983 <<mode.
mpAtom0->GetName()<<
"-"
4985 <<
"(Restraint="<<l-d0<<
"s"<<sigma<<
"d"<<delta<<
"):"
4986 <<l<<
"->"<<l+change<<endl;
4990 else change=(2.*(REAL)rand()-(REAL)RAND_MAX)/(REAL)RAND_MAX*amplitude*0.1;
4999 const bool respectRestraint)
5008 const REAL vx=dy10*dz12-dz10*dy12;
5009 const REAL vy=dz10*dx12-dx10*dz12;
5010 const REAL vz=dx10*dy12-dy10*dx12;
5013 if((abs(vx)+abs(vy)+abs(vz))<1e-6)
return change;
5017 const REAL norm= sqrt( (dx10*dx10+dy10*dy10+dz10*dz10)*(dx12*dx12+dy12*dy12+dz12*dz12)+1e-6);
5018 angle0=(dx10*dx12+dy10*dy12+dz10*dz12)/norm;
5019 if(angle0>=1) angle0=0;
5022 if(angle0<=-1) angle0=M_PI;
5023 else angle0= acos(angle0);
5026 const REAL a0=angle0-mode.
mpBondAngle->GetAngle0();
5027 const REAL sigma=mode.
mpBondAngle->GetAngleSigma();
5028 const REAL delta=mode.
mpBondAngle->GetAngleDelta();
5031 REAL a1=a0+(REAL)(2*rand()-RAND_MAX)/(REAL)RAND_MAX*amplitude*mode.
mBaseAmplitude;
5032 if(a1> delta)a1= delta;
5033 if(a1<-delta)a1=-delta;
5037 if((a0+change)>(delta+sigma*5.0)) change= delta+sigma*5.0-a0;
5038 else if((a0+change)<(-delta-sigma*5.0)) change=-delta-sigma*5.0-a0;
5042 cout<<
"ANGLE change("<<change*RAD2DEG<<
"):"
5043 <<mode.
mpAtom0->GetName()<<
"-"
5044 <<mode.
mpAtom1->GetName()<<
"-"
5046 <<
"(Restraint="<<(angle0-a0)*RAD2DEG<<
"s"<<sigma*RAD2DEG<<
"d"<<delta*RAD2DEG<<
"):"
5047 <<angle0*RAD2DEG<<
"->"<<(angle0+change)*RAD2DEG<<endl;
5051 else change=(2.*(REAL)rand()-(REAL)RAND_MAX)/(REAL)RAND_MAX*mode.
mBaseAmplitude*amplitude;
5056 const bool respectRestraint)
5062 if((abs(dx)+abs(dy)+abs(dz))<1e-6)
return change;
5071 REAL a1=a0+(REAL)(2*rand()-RAND_MAX)/(REAL)RAND_MAX*amplitude*mode.
mBaseAmplitude;
5072 if(a1> delta)a1= delta;
5073 if(a1<-delta)a1=-delta;
5077 if((a0+change)>(delta+sigma*5.0)) change= delta+sigma*5.0-a0;
5078 else if((a0+change)<(-delta-sigma*5.0)) change=-delta-sigma*5.0-a0;
5082 cout<<
"TORSION change ("
5083 <<mode.
mpAtom1->GetName()<<
"-"<<mode.
mpAtom2->GetName()<<
"):"<<endl
5084 <<
" initial angle="<<angle0*RAD2DEG
5086 <<(angle0-a0)*RAD2DEG<<
"s"<<sigma*RAD2DEG<<
"d"<<delta*RAD2DEG<<
"):"
5087 <<endl<<
" New angle:"<<(angle0+change)*RAD2DEG<<
", change="<<change*RAD2DEG<<endl;
5091 else change=(REAL)(2.*rand()-RAND_MAX)/(REAL)RAND_MAX*mode.
mBaseAmplitude*amplitude;
5104 mClockAtomPosition.
Click();
5108 void BuildZMatrixRecursive(
long &z,
const long curr,
5109 const vector<MolAtom*> &vpAtom,
5110 const map<
MolAtom *, set<MolAtom *> > &connT,
5111 vector<MolZAtom> &zmatrix,
5112 const map<const MolAtom*,long> &vIndex,
5113 vector<long> &vZIndex,
5114 vector<long> &vrZIndex)
5116 zmatrix[z].mpPow=&(vpAtom[curr]->GetScatteringPower());
5120 map<MolAtom *, set<MolAtom *> >::const_iterator pConn=connT.find(vpAtom[curr]);
5121 const long nc=pConn->second.size();
5122 vector<long> conn(nc);
5123 vector<long> zconn(nc);
5124 vector<long>::iterator pos=conn.begin();
5125 vector<long>::iterator zpos=zconn.begin();
5126 for(set<MolAtom *>::const_iterator pos1=pConn->second.begin();pos1!=pConn->second.end();++pos1)
5128 *pos = vIndex.find(*pos1)->second;
5129 *zpos = vZIndex[*pos];
5130 cout<<(*pos1)->GetName()<<
"("<<*pos<<
","<<*zpos<<
")"<<endl;
5133 sort(conn.begin(),conn.end());
5134 sort(zconn.begin(),zconn.end());
5138 const long b=zconn[nc-1];
5139 zmatrix[z].mBondAtom=b;
5140 zmatrix[z].mBondLength=
GetBondLength(*vpAtom[vrZIndex[b]],*vpAtom[curr]);
5143 const long a=zmatrix[b].mBondAtom;
5144 zmatrix[z].mBondAngleAtom=a;
5145 zmatrix[z].mBondAngle=
GetBondAngle(*vpAtom[vrZIndex[a]],*vpAtom[vrZIndex[b]],*vpAtom[curr]);
5148 const long d=zmatrix[b].mBondAngleAtom;
5149 zmatrix[z].mDihedralAtom=d;
5150 zmatrix[z].mDihedralAngle=fmod(
GetDihedralAngle(*vpAtom[vrZIndex[d]],*vpAtom[vrZIndex[a]],
5151 *vpAtom[vrZIndex[b]],*vpAtom[curr])+2*M_PI,
5156 zmatrix[z].mDihedralAtom=0;
5157 zmatrix[z].mDihedralAngle=0;
5162 zmatrix[z].mBondAngleAtom=0;
5163 zmatrix[z].mBondAngle=0;
5168 zmatrix[z].mBondAtom=0;
5169 zmatrix[z].mBondLength=0;
5173 for(pos=conn.begin();pos!=conn.end();++pos)
5177 if(vZIndex[*pos]==-1)
5178 BuildZMatrixRecursive(z,*pos,vpAtom,connT,zmatrix,vIndex,vZIndex,vrZIndex);
5188 map<const MolAtom*,long> vIndex;
5191 for(vector<MolAtom*>::const_iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
5197 for(
long i=0;i<n;++i)
5205 for(set<MolAtom *>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
5206 if((vIndex[*pos]<i)&&(vIndex[*pos]>b)) b=vIndex[*pos];
5213 const long a= (b==0)?1 :
mAsZMatrix[b].mBondAtom;
5224 for(set<MolAtom *>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
5225 if((vIndex[*pos]<i) && (vIndex[*pos]!=b) && (vIndex[*pos]>d)) d=vIndex[*pos];
5231 for(set<MolAtom *>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
5232 if((vIndex[*pos]<i) && (vIndex[*pos]!=a) && (vIndex[*pos]>d)) d=vIndex[*pos];
5238 for(set<MolAtom *>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
5239 if( (vIndex[*pos]<i) && (vIndex[*pos]!=a)
5240 &&(vIndex[*pos]!=b) && (vIndex[*pos]>d)) d=vIndex[*pos];
5244 for(
long j=0;j<i;++j)
5245 if((j!=a) &&(j!=b) && (j>d)) d=j;
5261 vector<long> vZIndex(n);
5263 vector<long> vrZIndex(n);
5264 for(
long i=0;i<n;++i)
5290 const map<
MolAtom *, set<MolAtom *> > &connect,
5291 list<MolAtom *> &atomlist,
5292 map<set<MolAtom *>,list<MolAtom *> > &ringlist)
5294 list<MolAtom *>::const_iterator f=find(atomlist.begin(),atomlist.end(),currentAtom);
5295 if(f!=atomlist.end())
5298 cout<<currentAtom->GetName()<<
" was already in the list : ring found !"<<endl;
5299 for(list<MolAtom *>::const_iterator atom=atomlist.begin();atom!=atomlist.end();++atom)
5300 cout<<(*atom)->GetName()<<
" ";
5303 set<MolAtom *> ring1;
5304 list<MolAtom *> ring2;
5305 for(list<MolAtom *>::const_iterator pos=f;pos!=atomlist.end();++pos)
5308 ring2.push_back(*pos);
5310 ringlist.insert(make_pair(ring1,ring2));
5314 atomlist.push_back(currentAtom);
5315 map<MolAtom *,set<MolAtom *> >::const_iterator c=connect.find(currentAtom);
5316 if(c != connect.end())
5318 set<MolAtom *>::const_iterator pos;
5319 for(pos=c->second.begin();pos!=c->second.end();++pos)
5321 if(*pos==previousAtom)
continue;
5325 atomlist.pop_back();
5332 if(mClockRingList>mClockConnectivityTable)
return;
5333 VFN_DEBUG_ENTRY(
"Molecule::BuildRingList()",7)
5334 for(vector<MolAtom*>::const_iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
5335 (*pos)->SetIsInRing(
false);
5336 list<MolAtom *> atomlist;
5338 map<set<MolAtom *>,list<MolAtom *> > ringlist;
5339 for(
unsigned long i=0;i<
mvpAtom.size();i++)
5344 for(map<set<MolAtom *>,list<MolAtom *> >::const_iterator pos0=ringlist.begin();pos0!=ringlist.end();pos0++)
5347 std::list<MolAtom*> *pList=&(
mvRing.back().GetAtomList());
5349 cout<<
"Found ring:";
5351 for(list<MolAtom *>::const_iterator atom=pos0->second.begin();atom!=pos0->second.end();++atom)
5353 pList->push_back(*atom);
5354 (*atom)->SetIsInRing(
true);
5356 cout<<(*atom)->GetName()<<
" ";
5364 cout<<
"Rings found :"<<ringlist.size()<<
", "<<
mvRing.size()<<
" unique."<<endl;
5365 mClockRingList.
Click();
5366 VFN_DEBUG_EXIT(
"Molecule::BuildRingList()",7)
5371 if( (mClockConnectivityTable>mClockBondList)
5372 &&(mClockConnectivityTable>mClockAtomList))
return;
5373 VFN_DEBUG_ENTRY(
"Molecule::BuildConnectivityTable()",5)
5374 TAU_PROFILE(
"Molecule::BuildConnectivityTable()",
"void ()",TAU_DEFAULT);
5379 for (
unsigned long i = 0; i <
mvpAtom.size(); ++i)
5383 for(
unsigned long i=0;i<
mvpBond.size();++i)
5391 map<MolAtom *,set<MolAtom *> >::const_iterator pos;
5394 cout<<
"Atom "<<pos->first->GetName()<<
" is connected to atoms: ";
5395 set<MolAtom *>::const_iterator pos1;
5396 for(pos1=pos->second.begin();pos1!=pos->second.end();++pos1)
5398 cout<<(*pos1)->GetName()<<
" ";
5401 if(pos->second.size()>24)
5402 throw ObjCrystException(
"Molecule: one atom ("+pos->first->GetName()+
") has more than 24 bonds !");
5408 for (
unsigned long i = 0; i <
mvpAtom.size(); ++i)
5410 cout <<
"Warning: Molecule '" << this->
GetName() <<
"' ConnectivityTable: Atom #"
5411 << i <<
"(" <<
mvpAtom[i]->GetName() <<
") has no bond !" << endl;
5413 mClockConnectivityTable.
Click();
5414 VFN_DEBUG_EXIT(
"Molecule::BuildConnectivityTable()",5)
5418 mpAtom1(&at1),mpAtom2(&at2),mBaseRotationAmplitude(M_PI*0.04)
5423 if( (mClockRotorGroup>mClockBondList)
5424 &&(mClockRotorGroup>mClockAtomList)
5425 &&(mClockRotorGroup>mClockBondAngleList)
5426 &&(mClockRotorGroup>mClockDihedralAngleList))
return;
5427 VFN_DEBUG_ENTRY(
"Molecule::BuildRotorGroup()",5)
5428 TAU_PROFILE(
"Molecule::BuildRotorGroup()",
"void ()",TAU_DEFAULT);
5435 for(
unsigned long i=0;i<
mvpBond.size();++i)
5439 *
const atom2=&(
mvpBond[i]->GetAtom2());
5440 for(
unsigned int j=1;j<=2;++j)
5442 const set<MolAtom*> *pConn;
5451 for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
5453 if((j==1)&&(*pos==atom2))
continue;
5454 if((j==2)&&(*pos==atom1))
continue;
5485 for(vector<MolAtom*>::const_iterator atom1=this->GetAtomList().begin();
5486 atom1!=this->GetAtomList().end();++atom1)
5489 vector<MolAtom*>::const_iterator atom2=atom1;
5491 for(;atom2!=this->GetAtomList().end();++atom2)
5493 for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
5495 if(*pos==*atom2)
continue;
5502 set<MolAtom*>::const_iterator check
5523 for(
unsigned int i=1;i<=3;++i)
5525 list<RotorGroup> *pRotorGroup1;
5532 for(list<RotorGroup>::iterator pos1=pRotorGroup1->begin();
5533 pos1!=pRotorGroup1->end();++pos1)
5535 for(
unsigned int j=i;j<=3;++j)
5537 list<RotorGroup> *pRotorGroup2;
5544 for(list<RotorGroup>::iterator pos2=pRotorGroup2->begin();
5545 pos2!=pRotorGroup2->end();)
5547 if(pos2==pos1) {++pos2;
continue;}
5548 if(( ((pos1->mpAtom1 == pos2->mpAtom1) && (pos1->mpAtom2 == pos2->mpAtom2))
5549 ||((pos1->mpAtom2 == pos2->mpAtom1) && (pos1->mpAtom1 == pos2->mpAtom2)))
5550 &&pos1->mvRotatedAtomList.size() == pos2->mvRotatedAtomList.size())
5553 for(set<MolAtom*>::const_iterator pos=pos1->mvRotatedAtomList.begin();
5554 pos!=pos1->mvRotatedAtomList.end();++pos)
5556 set<MolAtom*>::const_iterator tmp=pos2->mvRotatedAtomList.find(*pos);
5557 if(tmp == pos2->mvRotatedAtomList.end())
5566 cout<<
"Identical groups:"<<endl;
5568 <<pos1->mpAtom1->GetName()<<
"-"
5569 <<pos1->mpAtom2->GetName()<<
" : ";
5570 for(set<MolAtom*>::iterator pos=pos1->mvRotatedAtomList.begin();
5571 pos!=pos1->mvRotatedAtomList.end();++pos)
5572 cout<<(*pos)->GetName()<<
" ";
5575 <<pos2->mpAtom1->GetName()<<
"-"
5576 <<pos2->mpAtom2->GetName()<<
" : ";
5577 for(set<MolAtom*>::iterator pos=pos2->mvRotatedAtomList.begin();
5578 pos!=pos2->mvRotatedAtomList.end();++pos)
5579 cout<<(*pos)->GetName()<<
" ";
5582 pos2=pRotorGroup2->erase(pos2);
5594 for(
unsigned int i=1;i<=3;++i)
5596 list<RotorGroup> *pRotorGroup1;
5603 for(list<RotorGroup>::iterator pos=pRotorGroup1->begin();
5604 pos!=pRotorGroup1->end();)
5607 for(
unsigned int j=0;j<36;++j)
5609 const REAL angle=(REAL)j*M_PI/36.;
5611 pos->mvRotatedAtomList,angle);
5619 case 1: cout<<
"Rotation Group around bond :";
break;
5620 case 2: cout<<
"Rotation Group (single chain) around bond :";
break;
5621 case 3: cout<<
"Rotation Group (internal) between :";
break;
5623 cout <<pos->mpAtom1->GetName()<<
"-"
5624 <<pos->mpAtom2->GetName()<<
" : ";
5625 for(set<MolAtom *>::iterator pos1=pos->mvRotatedAtomList.begin();
5626 pos1!=pos->mvRotatedAtomList.end();++pos1)
5627 cout<<(*pos1)->GetName()<<
" ";
5628 cout<<
" <d(LLK)>="<< llk/36.;
5632 pos = pRotorGroup1->erase(pos);
5643 for(vector<MolBond*>::iterator pos=
mvpBond.begin();pos!=
mvpBond.end();++pos)
5644 (*pos)->SetFreeTorsion(
false);
5648 vector<MolBond*>::iterator bd=this->
FindBond((*pos->mpAtom1),(*pos->mpAtom2));
5649 if(bd!=
mvpBond.end()) (*bd)->SetFreeTorsion(
true);
5652 mClockRotorGroup.
Click();
5653 VFN_DEBUG_EXIT(
"Molecule::BuildRotorGroup()",5)
5659 if( (mClockStretchModeBondLength>mClockBondList)
5660 &&(mClockStretchModeBondLength>mClockAtomList)
5661 &&(mClockStretchModeBondLength>mClockBondAngleList)
5662 &&(mClockStretchModeBondLength>mClockDihedralAngleList))
return;
5664 VFN_DEBUG_ENTRY(
"Molecule::BuildStretchModeBondLength()",7)
5666 TAU_PROFILE(
"Molecule::BuildStretchModeBondLength()",
"void ()",TAU_DEFAULT);
5667 TAU_PROFILE_TIMER(timer1,
"Molecule::BuildStretchModeBondLength 1",
"", TAU_FIELD);
5668 TAU_PROFILE_TIMER(timer2,
"Molecule::BuildStretchModeBondLength 2",
"", TAU_FIELD);
5669 TAU_PROFILE_TIMER(timer3,
"Molecule::BuildStretchModeBondLength 3",
"", TAU_FIELD);
5670 TAU_PROFILE_TIMER(timer4,
"Molecule::BuildStretchModeBondLength 4",
"", TAU_FIELD);
5671 TAU_PROFILE_TIMER(timer5,
"Molecule::BuildStretchModeBondLength 5",
"", TAU_FIELD);
5675 TAU_PROFILE_START(timer1);
5676 for(
unsigned long i=0;i<
mvpBond.size();++i)
5681 for(
unsigned int j=1;j<=2;++j)
5683 const set<MolAtom*> *pConn;
5698 for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
5700 if((j==1)&&(*pos==atom2))
continue;
5701 if((j==2)&&(*pos==atom1))
continue;
5707 if( ((j==1)&&(ct2>0)) || ((j==2)&&(ct1>0)) )
5726 TAU_PROFILE_STOP(timer1);
5732 TAU_PROFILE_START(timer5);
5734 for(vector<RigidGroup*>::const_iterator group=
mvRigidGroup.begin();
5738 for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
5739 ct += pos->mvTranslatedAtomList.count(*at);
5740 if((ct>0)&&(ct!=(*group)->size()))
5744 cout<<
" Breaks Rigid Group:";
5745 for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
5746 cout<<(*at)->GetName()<<
" ";
5754 TAU_PROFILE_STOP(timer5);
5758 unsigned long paramSetRandom[5];
5759 for(
unsigned long i=0;i<5;++i)
5761 for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
5763 (*pos)->SetX(100.*rand()/(REAL) RAND_MAX);
5764 (*pos)->SetY(100.*rand()/(REAL) RAND_MAX);
5765 (*pos)->SetZ(100.*rand()/(REAL) RAND_MAX);
5773 TAU_PROFILE_START(timer2);
5775 pos->mvpBrokenBond.clear();
5776 for(vector<MolBond*>::const_iterator r=
mvpBond.begin();r!=
mvpBond.end();++r)
5779 for(set<MolAtom *>::const_iterator at=pos->mvTranslatedAtomList.begin();
5780 at!=pos->mvTranslatedAtomList.end();++at)
5782 if(*at==&((*r)->GetAtom1())) ct++;
5783 if(*at==&((*r)->GetAtom2())) ct++;
5786 if((ct!=0)&&(ct !=2)) pos->mvpBrokenBond.insert(make_pair(*r,0));
5790 int nb=pos->mvpBrokenBond.size();
5791 if(pos->mpBond!=0) nb -= 1;
5792 if(nb>0) keep=
false;
5796 TAU_PROFILE_STOP(timer2);
5802 TAU_PROFILE_START(timer3);
5804 pos->mvpBrokenBondAngle.clear();
5808 for(set<MolAtom *>::const_iterator at=pos->mvTranslatedAtomList.begin();
5809 at!=pos->mvTranslatedAtomList.end();++at)
5811 if(*at==&((*r)->GetAtom1())) ct++;
5812 if(*at==&((*r)->GetAtom2())) ct++;
5813 if(*at==&((*r)->GetAtom3())) ct++;
5816 if((ct==0)||(ct==3)) broken=
false;
5820 for(
unsigned long i=0;i<5;++i)
5823 pos->CalcDeriv(
false);
5824 (*r)->GetLogLikelihood(
true,
true);
5825 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
5828 if(abs(d)<=0.01) broken=
false;
5830 if(broken) pos->mvpBrokenBondAngle.insert(make_pair(*r,0));
5834 if(pos->mvpBrokenBondAngle.size()>0) keep=
false;
5838 TAU_PROFILE_STOP(timer3);
5844 TAU_PROFILE_START(timer4);
5846 pos->mvpBrokenDihedralAngle.clear();
5850 for(set<MolAtom *>::const_iterator at=pos->mvTranslatedAtomList.begin();
5851 at!=pos->mvTranslatedAtomList.end();++at)
5853 if(*at==&((*r)->GetAtom1())) ct++;
5854 if(*at==&((*r)->GetAtom2())) ct++;
5855 if(*at==&((*r)->GetAtom3())) ct++;
5856 if(*at==&((*r)->GetAtom4())) ct++;
5859 if((ct==0)||(ct==4)) broken=
false;
5863 for(
unsigned long i=0;i<5;++i)
5866 pos->CalcDeriv(
false);
5867 (*r)->GetLogLikelihood(
true,
true);
5868 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
5871 if(abs(d)<=0.01) broken=
false;
5873 if(broken) pos->mvpBrokenDihedralAngle.insert(make_pair(*r,0));
5877 if(pos->mvpBrokenDihedralAngle.size()>0) keep=
false;
5881 TAU_PROFILE_STOP(timer4);
5884 for(
unsigned long i=0;i<5;++i) this->
ClearParamSet(paramSetRandom[i]);
5886 cout<<
"List of Bond Length stretch modes"<<endl;
5890 cout<<
" Bond:"<<pos->mpAtom0->GetName()<<
"-"<<pos->mpAtom1->GetName()<<
", Translated Atoms: ";
5891 for(set<MolAtom*>::const_iterator atom=pos->mvTranslatedAtomList.begin();
5892 atom!=pos->mvTranslatedAtomList.end();++atom)
5894 cout<<(*atom)->GetName()<<
",";
5896 if(pos->mpBond!=0) cout<<
" ; restrained to length="<<pos->mpBond->GetLength0()
5897 <<
", sigma="<<pos->mpBond->GetLengthSigma()
5898 <<
", delta="<<pos->mpBond->GetLengthDelta();
5899 if(pos->mvpBrokenBond.size()>0)
5901 cout<<endl<<
" Broken bonds:";
5902 for(map<const MolBond*,REAL>::const_iterator bond=pos->mvpBrokenBond.begin();
5903 bond!=pos->mvpBrokenBond.end();++bond)
5904 cout<<bond->first->GetName()<<
", ";
5906 if(pos->mvpBrokenBondAngle.size()>0)
5908 cout<<endl<<
" Broken bond angles:";
5909 for(map<const MolBondAngle*,REAL>::const_iterator angle=pos->mvpBrokenBondAngle.begin();
5910 angle!=pos->mvpBrokenBondAngle.end();++angle)
5911 cout<<angle->first->GetName()<<
", ";
5913 if(pos->mvpBrokenDihedralAngle.size()>0)
5915 cout<<endl<<
" Broken dihedral angles:";
5916 for(map<const MolDihedralAngle*,REAL>::const_iterator
5917 angle=pos->mvpBrokenDihedralAngle.begin();
5918 angle!=pos->mvpBrokenDihedralAngle.end();++angle)
5919 cout<<angle->first->GetName()<<
", ";
5924 mClockStretchModeBondLength.
Click();
5925 VFN_DEBUG_EXIT(
"Molecule::BuildStretchModeBondLength()",7)
5931 if( (mClockStretchModeBondAngle>mClockBondList)
5932 &&(mClockStretchModeBondAngle>mClockAtomList)
5933 &&(mClockStretchModeBondAngle>mClockBondAngleList)
5934 &&(mClockStretchModeBondAngle>mClockDihedralAngleList))
return;
5936 VFN_DEBUG_ENTRY(
"Molecule::BuildStretchModeBondAngle()",10)
5938 TAU_PROFILE(
"Molecule::BuildStretchModeBondAngle()",
"void ()",TAU_DEFAULT);
5939 TAU_PROFILE_TIMER(timer1,
"Molecule::BuildStretchModeBondAngle 1",
"", TAU_FIELD);
5940 TAU_PROFILE_TIMER(timer2,
"Molecule::BuildStretchModeBondAngle 2",
"", TAU_FIELD);
5941 TAU_PROFILE_TIMER(timer3,
"Molecule::BuildStretchModeBondAngle 3",
"", TAU_FIELD);
5942 TAU_PROFILE_TIMER(timer4,
"Molecule::BuildStretchModeBondAngle 4",
"", TAU_FIELD);
5943 TAU_PROFILE_TIMER(timer5,
"Molecule::BuildStretchModeBondAngle 5",
"", TAU_FIELD);
5947 TAU_PROFILE_START(timer1);
5948 for(
unsigned long i=0;i<
mvpAtom.size();++i)
5952 if(pConn0->size()<2)
continue;
5953 for(set<MolAtom*>::const_iterator pos1=pConn0->begin();pos1!=pConn0->end();++pos1)
5955 set<MolAtom*>::const_iterator pos2=pos1;
5957 for(;pos2!=pConn0->end();++pos2)
5959 VFN_DEBUG_MESSAGE(
"Molecule::BuildStretchModeBondAngle():"<<i<<
","<<*pos1<<
","<<*pos2,10)
5964 if(&((*pos)->GetAtom2())==
mvpAtom[i])
5966 if( ((&((*pos)->GetAtom1())==*pos1)&&(&((*pos)->GetAtom3())==*pos2))
5967 ||((&((*pos)->GetAtom1())==*pos2)&&(&((*pos)->GetAtom3())==*pos1)))
5974 for(
unsigned int j=1;j<=2;++j)
5976 const set<MolAtom*> *pConn;
5997 for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
5999 if(*pos==
mvpAtom[i])
continue;
6039 TAU_PROFILE_STOP(timer1);
6044 TAU_PROFILE_START(timer5);
6046 for(vector<RigidGroup*>::const_iterator group=
mvRigidGroup.begin();
6050 for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
6051 ct += pos->mvRotatedAtomList.count(*at);
6055 ct += (*group)->count(pos->mpAtom1);
6056 if(ct!=(*group)->size())
6061 cout<<
" Breaks Rigid Group:";
6062 for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
6063 cout<<(*at)->GetName()<<
" ";
6072 TAU_PROFILE_STOP(timer5);
6076 unsigned long paramSetRandom[5];
6077 for(
unsigned long i=0;i<5;++i)
6079 for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
6081 (*pos)->SetX(100.*rand()/(REAL) RAND_MAX);
6082 (*pos)->SetY(100.*rand()/(REAL) RAND_MAX);
6083 (*pos)->SetZ(100.*rand()/(REAL) RAND_MAX);
6091 TAU_PROFILE_START(timer2);
6093 pos->mvpBrokenBond.clear();
6094 for(vector<MolBond*>::const_iterator r=
mvpBond.begin();r!=
mvpBond.end();++r)
6097 for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
6098 at!=pos->mvRotatedAtomList.end();++at)
6100 if(*at==&((*r)->GetAtom1())) ct++;
6101 if(*at==&((*r)->GetAtom2())) ct++;
6105 if((ct==0)||(ct==2)) broken=
false;
6109 for(
unsigned long i=0;i<5;++i)
6112 pos->CalcDeriv(
false);
6113 (*r)->GetLogLikelihood(
true,
true);
6114 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6117 if(abs(d)<=0.01) broken=
false;
6119 if(broken) pos->mvpBrokenBond.insert(make_pair(*r,0));
6123 if(pos->mvpBrokenBond.size()>0) keep=
false;
6127 TAU_PROFILE_STOP(timer2);
6133 TAU_PROFILE_START(timer3);
6135 pos->mvpBrokenBondAngle.clear();
6139 for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
6140 at!=pos->mvRotatedAtomList.end();++at)
6142 if(*at==&((*r)->GetAtom1())) ct++;
6143 if(*at==&((*r)->GetAtom2())) ct++;
6144 if(*at==&((*r)->GetAtom3())) ct++;
6147 if(ct==0) broken=
false;
6148 if(ct==3) broken=
false;
6152 for(
unsigned long i=0;i<5;++i)
6155 pos->CalcDeriv(
false);
6156 (*r)->GetLogLikelihood(
true,
true);
6157 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6160 if(abs(d)<=0.01) broken=
false;
6162 if(broken) pos->mvpBrokenBondAngle.insert(make_pair(*r,0));
6166 int nb=pos->mvpBrokenBond.size();
6167 if(pos->mpBondAngle!=0) nb -= 1;
6168 if(nb>0) keep=
false;
6172 TAU_PROFILE_STOP(timer3);
6178 TAU_PROFILE_START(timer4);
6180 pos->mvpBrokenDihedralAngle.clear();
6184 for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
6185 at!=pos->mvRotatedAtomList.end();++at)
6187 if(*at==&((*r)->GetAtom1())) ct++;
6188 if(*at==&((*r)->GetAtom2())) ct++;
6189 if(*at==&((*r)->GetAtom3())) ct++;
6190 if(*at==&((*r)->GetAtom4())) ct++;
6193 if(ct==0) broken=
false;
6194 if(ct==4) broken=
false;
6198 for(
unsigned long i=0;i<5;++i)
6201 pos->CalcDeriv(
false);
6202 (*r)->GetLogLikelihood(
true,
true);
6203 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6206 if(abs(d)<=0.01) broken=
false;
6208 if(broken) pos->mvpBrokenDihedralAngle.insert(make_pair(*r,0));
6212 if(pos->mvpBrokenDihedralAngle.size()>0) keep=
false;
6216 TAU_PROFILE_STOP(timer4);
6219 for(
unsigned long i=0;i<5;++i) this->
ClearParamSet(paramSetRandom[i]);
6221 cout<<
"List of Bond Angle stretch modes"<<endl;
6225 cout<<
" Angle:"<<pos->mpAtom0->GetName()<<
"-"
6226 <<pos->mpAtom1->GetName()<<
"-"
6227 <<pos->mpAtom2->GetName()<<
", Rotated Atoms: ";
6228 for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
6229 atom!=pos->mvRotatedAtomList.end();++atom)
6231 cout<<(*atom)->GetName()<<
",";
6233 if(pos->mpBondAngle!=0) cout<<
" ; restrained to angle="<<pos->mpBondAngle->GetAngle0()*RAD2DEG
6234 <<
"�, sigma="<<pos->mpBondAngle->GetAngleSigma()*RAD2DEG
6235 <<
"�, delta="<<pos->mpBondAngle->GetAngleDelta()*RAD2DEG<<
"�";
6236 if(pos->mvpBrokenBond.size()>0)
6238 cout<<endl<<
" Broken bonds:";
6239 for(map<const MolBond*,REAL>::const_iterator bond=pos->mvpBrokenBond.begin();
6240 bond!=pos->mvpBrokenBond.end();++bond)
6241 cout<<bond->first->GetName()<<
", ";
6243 if(pos->mvpBrokenBondAngle.size()>0)
6245 cout<<endl<<
" Broken bond angles:";
6246 for(map<const MolBondAngle*,REAL>::const_iterator angle=pos->mvpBrokenBondAngle.begin();
6247 angle!=pos->mvpBrokenBondAngle.end();++angle)
6248 cout<<angle->first->GetName()<<
", ";
6250 if(pos->mvpBrokenDihedralAngle.size()>0)
6252 cout<<endl<<
" Broken dihedral angles:";
6253 for(map<const MolDihedralAngle*,REAL>::const_iterator
6254 angle=pos->mvpBrokenDihedralAngle.begin();
6255 angle!=pos->mvpBrokenDihedralAngle.end();++angle)
6256 cout<<angle->first->GetName()<<
", ";
6261 mClockStretchModeBondAngle.
Click();
6262 VFN_DEBUG_EXIT(
"Molecule::BuildStretchModeBondAngle()",10)
6268 if( (mClockStretchModeTorsion>mClockBondList)
6269 &&(mClockStretchModeTorsion>mClockAtomList)
6270 &&(mClockStretchModeTorsion>mClockBondAngleList)
6271 &&(mClockStretchModeTorsion>mClockDihedralAngleList))
return;
6273 VFN_DEBUG_ENTRY(
"Molecule::BuildStretchModeTorsion()",7)
6274 TAU_PROFILE(
"Molecule::BuildStretchModeTorsion()",
"void ()",TAU_DEFAULT);
6275 TAU_PROFILE_TIMER(timer1,
"Molecule::BuildStretchModeTorsion 1",
"", TAU_FIELD);
6276 TAU_PROFILE_TIMER(timer2,
"Molecule::BuildStretchModeTorsion 2",
"", TAU_FIELD);
6277 TAU_PROFILE_TIMER(timer3,
"Molecule::BuildStretchModeTorsion 3",
"", TAU_FIELD);
6278 TAU_PROFILE_TIMER(timer4,
"Molecule::BuildStretchModeTorsion 4",
"", TAU_FIELD);
6279 TAU_PROFILE_TIMER(timer5,
"Molecule::BuildStretchModeTorsion 5",
"", TAU_FIELD);
6284 for(
unsigned long i=0;i<
mvpBond.size();++i)
6289 for(
unsigned int j=1;j<=2;++j)
6291 const set<MolAtom*> *pConn;
6298 for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
6300 if((*pos==atom2)||(*pos==atom1))
continue;
6312 if( ((&((*dih)->GetAtom2())==atom1) && (&((*dih)->GetAtom3())==atom2))
6313 ||((&((*dih)->GetAtom3())==atom1) && (&((*dih)->GetAtom2())==atom2)))
6315 const unsigned long ct1=
mvStretchModeTorsion.back().mvRotatedAtomList.count(&((*dih)->GetAtom1()));
6316 const unsigned long ct4=
mvStretchModeTorsion.back().mvRotatedAtomList.count(&((*dih)->GetAtom4()));
6344 if( ( ((mode->mpAtom1==atom1)&&(mode->mpAtom2==atom2))
6345 ||((mode->mpAtom1==atom2)&&(mode->mpAtom2==atom1)))
6367 for(
unsigned long i=0;i<
mvpBond.size();++i)
6372 for(
unsigned int j=1;j<=2;++j)
6374 const set<MolAtom*> *pConn;
6378 for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
6380 if((*pos==atom2)||(*pos==atom1))
continue;
6394 if( ((&((*dih)->GetAtom2())==atom1) && (&((*dih)->GetAtom3())==atom2))
6395 ||((&((*dih)->GetAtom3())==atom1) && (&((*dih)->GetAtom2())==atom2)))
6397 const unsigned long ct1=
mvStretchModeTorsion.back().mvRotatedAtomList.count(&((*dih)->GetAtom1()));
6398 const unsigned long ct4=
mvStretchModeTorsion.back().mvRotatedAtomList.count(&((*dih)->GetAtom4()));
6426 if( ( ((mode->mpAtom1==atom1)&&(mode->mpAtom2==atom2))
6427 ||((mode->mpAtom1==atom2)&&(mode->mpAtom2==atom1)))
6455 TAU_PROFILE_START(timer5);
6457 for(vector<RigidGroup*>::const_iterator group=
mvRigidGroup.begin();
6461 for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
6462 ct += pos->mvRotatedAtomList.count(*at);
6466 ct += (*group)->count(pos->mpAtom1);
6467 ct += (*group)->count(pos->mpAtom2);
6468 if(ct!=(*group)->size())
6473 cout<<
" Breaks Rigid Group:";
6474 for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
6475 cout<<(*at)->GetName()<<
" ";
6484 TAU_PROFILE_STOP(timer5);
6488 unsigned long paramSetRandom[5];
6489 for(
unsigned long i=0;i<5;++i)
6491 for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
6493 (*pos)->SetX(100.*rand()/(REAL) RAND_MAX);
6494 (*pos)->SetY(100.*rand()/(REAL) RAND_MAX);
6495 (*pos)->SetZ(100.*rand()/(REAL) RAND_MAX);
6503 TAU_PROFILE_START(timer2);
6505 pos->mvpBrokenBond.clear();
6506 for(vector<MolBond*>::const_iterator r=
mvpBond.begin();r!=
mvpBond.end();++r)
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++;
6517 if((ct==0)||(ct==2)) broken=
false;
6521 for(
unsigned long i=0;i<5;++i)
6524 pos->CalcDeriv(
false);
6525 (*r)->GetLogLikelihood(
true,
true);
6526 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6529 if(abs(d)<=0.01) broken=
false;
6531 if(broken) pos->mvpBrokenBond.insert(make_pair(*r,0));
6535 if(pos->mvpBrokenBond.size()>0) keep=
false;
6539 TAU_PROFILE_STOP(timer2);
6545 TAU_PROFILE_START(timer3);
6547 pos->mvpBrokenBondAngle.clear();
6551 for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
6552 at!=pos->mvRotatedAtomList.end();++at)
6554 if(*at==&((*r)->GetAtom1())) ct++;
6555 if(*at==&((*r)->GetAtom2())) ct++;
6556 if(*at==&((*r)->GetAtom3())) ct++;
6559 if((ct==0)||(ct==3)) broken=
false;
6563 for(
unsigned long i=0;i<5;++i)
6566 pos->CalcDeriv(
false);
6567 (*r)->GetLogLikelihood(
true,
true);
6568 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6571 if(abs(d)<=0.01) broken=
false;
6573 if(broken) pos->mvpBrokenBondAngle.insert(make_pair(*r,0));
6577 if(pos->mvpBrokenBond.size()>0) keep=
false;
6581 TAU_PROFILE_STOP(timer3);
6587 TAU_PROFILE_START(timer4);
6589 pos->mvpBrokenDihedralAngle.clear();
6593 for(set<MolAtom *>::const_iterator at=pos->mvRotatedAtomList.begin();
6594 at!=pos->mvRotatedAtomList.end();++at)
6596 if(*at==&((*r)->GetAtom1())) ct++;
6597 if(*at==&((*r)->GetAtom2())) ct++;
6598 if(*at==&((*r)->GetAtom3())) ct++;
6599 if(*at==&((*r)->GetAtom4())) ct++;
6602 if((ct==0)||(ct==4)) broken=
false;
6606 for(
unsigned long i=0;i<5;++i)
6609 pos->CalcDeriv(
false);
6610 (*r)->GetLogLikelihood(
true,
true);
6611 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6614 if(abs(d)<=0.01) broken=
false;
6616 if(broken) pos->mvpBrokenDihedralAngle.insert(make_pair(*r,0));
6620 int nb=pos->mvpBrokenDihedralAngle.size();
6621 if(pos->mpDihedralAngle!=0) nb -= 1;
6622 if(nb>0) keep=
false;
6626 TAU_PROFILE_STOP(timer4);
6629 for(
unsigned long i=0;i<5;++i) this->
ClearParamSet(paramSetRandom[i]);
6637 cout<<
" Dihedral Angle:"
6638 <<pos->mpAtom1->GetName()<<
"-"
6639 <<pos->mpAtom2->GetName()<<
", Rotated Atoms: ";
6640 for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
6641 atom!=pos->mvRotatedAtomList.end();++atom)
6643 cout<<(*atom)->GetName()<<
",";
6645 if(pos->mpDihedralAngle!=0)
6646 cout<<endl<<
" ->restrained by dihedral angle "<<pos->mpDihedralAngle->GetName()
6647 <<
"to :"<<pos->mpDihedralAngle->GetAngle0()*RAD2DEG
6648 <<
"�, sigma="<<pos->mpDihedralAngle->GetAngleSigma()*RAD2DEG
6649 <<
"�, delta="<<pos->mpDihedralAngle->GetAngleDelta()*RAD2DEG<<
"�";
6650 if(pos->mvpBrokenBond.size()>0)
6652 cout<<endl<<
" Broken bonds:";
6653 for(map<const MolBond*,REAL>::const_iterator bond=pos->mvpBrokenBond.begin();
6654 bond!=pos->mvpBrokenBond.end();++bond)
6655 cout<<bond->first->GetName()<<
", ";
6657 if(pos->mvpBrokenBondAngle.size()>0)
6659 cout<<endl<<
" Broken bond angles:";
6660 for(map<const MolBondAngle*,REAL>::const_iterator angle=pos->mvpBrokenBondAngle.begin();
6661 angle!=pos->mvpBrokenBondAngle.end();++angle)
6662 cout<<angle->first->GetName()<<
", ";
6664 if(pos->mvpBrokenDihedralAngle.size()>0)
6666 cout<<endl<<
" Broken dihedral angles:";
6667 for(map<const MolDihedralAngle*,REAL>::const_iterator
6668 angle=pos->mvpBrokenDihedralAngle.begin();
6669 angle!=pos->mvpBrokenDihedralAngle.end();++angle)
6670 cout<<angle->first->GetName()<<
", ";
6675 mClockStretchModeTorsion.
Click();
6676 VFN_DEBUG_EXIT(
"Molecule::BuildStretchModeTorsion()",7)
6682 if( (mClockStretchModeTwist>mClockBondList)
6683 &&(mClockStretchModeTwist>mClockAtomList)
6684 &&(mClockStretchModeTwist>mClockBondAngleList)
6685 &&(mClockStretchModeTwist>mClockDihedralAngleList))
return;
6687 VFN_DEBUG_ENTRY(
"Molecule::BuildStretchModeTwist()",7)
6692 for(vector<MolAtom*>::const_iterator atom1=this->GetAtomList().begin();
6693 atom1!=this->GetAtomList().end();++atom1)
6696 vector<MolAtom*>::const_iterator atom2=atom1;
6698 for(;atom2!=this->GetAtomList().end();++atom2)
6700 for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
6702 if(*pos==*atom2)
continue;
6709 set<MolAtom*>::const_iterator check
6727 cout<<
"Rejecting StretchModeTwist ";
mvStretchModeTwist.back().Print(cout);cout<<endl;
6742 if( ( ((mode->mpAtom1==*atom1)&&(mode->mpAtom2==*atom2))
6743 ||((mode->mpAtom1==*atom2)&&(mode->mpAtom2==*atom1)))
6747 cout<<
"Duplicate StretchModeTwist ";
mvStretchModeTwist.back().Print(cout);cout<<endl;
6761 list<StretchModeTorsion>::reverse_iterator mode;
6764 if( ( ((mode->mpAtom1==*atom1)&&(mode->mpAtom2==*atom2))
6765 ||((mode->mpAtom1==*atom2)&&(mode->mpAtom2==*atom1)))
6769 cout<<
"Duplicate StretchModeTwist (with Torsion) ";
mvStretchModeTwist.back().Print(cout);cout<<endl;
6784 unsigned long paramSetRandom[5];
6785 for(
unsigned long i=0;i<5;++i)
6787 for(vector<MolAtom*>::iterator pos=
mvpAtom.begin();pos!=
mvpAtom.end();++pos)
6789 (*pos)->SetX(100.*rand()/(REAL) RAND_MAX);
6790 (*pos)->SetY(100.*rand()/(REAL) RAND_MAX);
6791 (*pos)->SetZ(100.*rand()/(REAL) RAND_MAX);
6799 pos->mvpBrokenBond.clear();
6800 pos->mvpBrokenBondAngle.clear();
6801 pos->mvpBrokenDihedralAngle.clear();
6804 cout<<
" DerivLLK for Twist mode around:"
6805 <<pos->mpAtom1->GetName()<<
"-"
6806 <<pos->mpAtom2->GetName()<<
": Moving atoms:";
6807 for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
6808 atom!=pos->mvRotatedAtomList.end();++atom)
6809 cout<<(*atom)->GetName()<<
",";
6812 for(vector<MolBond*>::const_iterator r=
mvpBond.begin();r!=
mvpBond.end();++r)
6814 (*r)->GetLogLikelihood(
true,
true);
6816 for(
unsigned long i=0;i<5;++i)
6819 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6824 cout<<
" Bond "<<(*r)->GetName()
6825 <<
": dLength/da="<<d<<endl;
6827 pos->mvpBrokenBond.insert(make_pair(*r,0.0));
6832 (*r)->GetLogLikelihood(
true,
true);
6834 for(
unsigned long i=0;i<5;++i)
6837 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6842 cout<<
" BondAngle:"<<(*r)->GetName()<<
": dAngle/da="<<d<<endl;
6844 pos->mvpBrokenBondAngle.insert(make_pair(*r,0.0));
6851 (*r)->GetLogLikelihood(
true,
true);
6853 for(
unsigned long i=0;i<5;++i)
6856 d += abs((*r)->GetDeriv(pos->mDerivXYZ));
6861 cout<<
" DihedralAngle:"<<(*r)->GetName()<<
": dAngle/da="<<d<<endl;
6863 pos->mvpBrokenDihedralAngle.insert(make_pair(*r,0.0));
6868 for(vector<RigidGroup*>::const_iterator group=
mvRigidGroup.begin();
6872 for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
6873 ct += pos->mvRotatedAtomList.count(*at);
6878 ct += (*group)->count(pos->mpAtom1);
6879 ct += (*group)->count(pos->mpAtom2);
6880 if(ct!=(*group)->size())
6884 cout<<
" Breaks Rigid Group:"<<ct<<
"!="<<(*group)->size()<<
":";
6885 for(set<MolAtom *>::const_iterator at=(*group)->begin();at!=(*group)->end();++at)
6886 cout<<(*at)->GetName()<<
" ";
6895 if(pos->mvpBrokenBond.size()+pos->mvpBrokenBondAngle.size()+pos->mvpBrokenDihedralAngle.size()) keep=
false;
6900 for(
unsigned long i=0;i<5;++i) this->
ClearParamSet(paramSetRandom[i]);
6908 cout<<
" Twist mode:"
6909 <<pos->mpAtom1->GetName()<<
"-"
6910 <<pos->mpAtom2->GetName()<<
", Rotated Atoms: ";
6911 for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
6912 atom!=pos->mvRotatedAtomList.end();++atom)
6914 cout<<(*atom)->GetName()<<
",";
6916 if(pos->mvpBrokenBond.size()>0)
6918 cout<<endl<<
" Broken bonds:";
6919 for(map<const MolBond*,REAL>::const_iterator bond=pos->mvpBrokenBond.begin();
6920 bond!=pos->mvpBrokenBond.end();++bond)
6921 cout<<bond->first->GetName()<<
", ";
6923 if(pos->mvpBrokenBondAngle.size()>0)
6925 cout<<endl<<
" Broken bond angles:";
6926 for(map<const MolBondAngle*,REAL>::const_iterator angle=pos->mvpBrokenBondAngle.begin();
6927 angle!=pos->mvpBrokenBondAngle.end();++angle)
6928 cout<<angle->first->GetName()<<
", ";
6930 if(pos->mvpBrokenDihedralAngle.size()>0)
6932 cout<<endl<<
" Broken dihedral angles:";
6933 for(map<const MolDihedralAngle*,REAL>::const_iterator
6934 angle=pos->mvpBrokenDihedralAngle.begin();
6935 angle!=pos->mvpBrokenDihedralAngle.end();++angle)
6936 cout<<angle->first->GetName()<<
", ";
6941 mClockStretchModeTwist.
Click();
6942 VFN_DEBUG_EXIT(
"Molecule::BuildStretchModeTwist()",7)
6947 VFN_DEBUG_ENTRY(
"Molecule::TuneGlobalOptimRotationAmplitude()",5)
6948 unsigned long initialConfig=this->
CreateParamSet(
"Initial Configuration");
6949 const unsigned int nbTest=100;
6958 REAL displacement=0;
6960 for(
unsigned int j=0;j<nbTest;j++)
6968 REAL xc=0.,yc=0.,zc=0.;
6971 x0[i]=
mvpAtom[i]->GetX(); xc += x0[i];
6972 y0[i]=
mvpAtom[i]->GetY(); yc += y0[i];
6973 z0[i]=
mvpAtom[i]->GetZ(); zc += z0[i];
6983 const REAL dx10=pos->mpAtom0->GetX()-pos->mpAtom1->GetX();
6984 const REAL dy10=pos->mpAtom0->GetY()-pos->mpAtom1->GetY();
6985 const REAL dz10=pos->mpAtom0->GetZ()-pos->mpAtom1->GetZ();
6986 const REAL dx12=pos->mpAtom2->GetX()-pos->mpAtom1->GetX();
6987 const REAL dy12=pos->mpAtom2->GetY()-pos->mpAtom1->GetY();
6988 const REAL dz12=pos->mpAtom2->GetZ()-pos->mpAtom1->GetZ();
6990 const REAL vx=dy10*dz12-dz10*dy12;
6991 const REAL vy=dz10*dx12-dx10*dz12;
6992 const REAL vz=dx10*dy12-dy10*dx12;
6993 this->
RotateAtomGroup(*(pos->mpAtom1),vx,vy,vz,pos->mvRotatedAtomList,0.01,
false);
6999 pos->mBaseAmplitude+=sqrt(abs(dx*dx+dy*dy+dz*dz));
7001 this->
RotateAtomGroup(*(pos->mpAtom1),vx,vy,vz,pos->mvRotatedAtomList,-0.01,
false);
7006 this->
RotateAtomGroup(*(pos->mpAtom1),*(pos->mpAtom2),pos->mvRotatedAtomList,0.01,
false);
7012 pos->mBaseAmplitude+=sqrt(abs(dx*dx+dy*dy+dz*dz));
7014 this->
RotateAtomGroup(*(pos->mpAtom1),*(pos->mpAtom2),pos->mvRotatedAtomList,-0.01,
false);
7017 for(
unsigned int k=0;k<10;++k)
7020 (mBaseRotationAmplitude,(REAL)rand(),(REAL)rand(),(REAL)rand());
7030 displacement+=sqrt(abs(dx*dx+dy*dy+dz*dz));
7039 pos->mBaseAmplitude=0.1*0.01/pos->mBaseAmplitude;
7040 if(pos->mBaseAmplitude>.17) pos->mBaseAmplitude=.17;
7042 if((pos->mvpBrokenBond.size()+pos->mvpBrokenBondAngle.size()+pos->mvpBrokenDihedralAngle.size())>0)
7044 pos->mBaseAmplitude/=10;
7049 <<pos->mpAtom0->GetName()<<
"-"
7050 <<pos->mpAtom1->GetName()<<
"-"
7051 <<pos->mpAtom2->GetName()<<
":";
7052 for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
7053 atom!=pos->mvRotatedAtomList.end();++atom) cout<<(*atom)->GetName()<<
",";
7054 cout <<
": base rotation="<<pos->mBaseAmplitude*RAD2DEG<<
" free="<<free<<endl;
7062 pos->mBaseAmplitude=0.1*0.01/pos->mBaseAmplitude;
7063 if(pos->mBaseAmplitude>.17) pos->mBaseAmplitude=.17;
7065 if((pos->mvpBrokenBond.size()+pos->mvpBrokenBondAngle.size()+pos->mvpBrokenDihedralAngle.size())>0)
7067 pos->mBaseAmplitude/=10;
7072 <<pos->mpAtom1->GetName()<<
"-"
7073 <<pos->mpAtom2->GetName()<<
":";
7074 for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
7075 atom!=pos->mvRotatedAtomList.end();++atom) cout<<(*atom)->GetName()<<
",";
7076 cout <<
": base rotation="<<pos->mBaseAmplitude*RAD2DEG<<
" free="<<free<<endl;
7084 pos->mBaseAmplitude=0.1*0.01/pos->mBaseAmplitude;
7085 if(pos->mBaseAmplitude>.17) pos->mBaseAmplitude=.17;
7087 if((pos->mvpBrokenBond.size()+pos->mvpBrokenBondAngle.size()+pos->mvpBrokenDihedralAngle.size())>0)
7089 pos->mBaseAmplitude/=10;
7094 <<pos->mpAtom1->GetName()<<
"-"
7095 <<pos->mpAtom2->GetName()<<
":";
7096 for(set<MolAtom*>::const_iterator atom=pos->mvRotatedAtomList.begin();
7097 atom!=pos->mvRotatedAtomList.end();++atom) cout<<(*atom)->GetName()<<
",";
7098 cout <<
": base rotation="<<pos->mBaseAmplitude*RAD2DEG<<
" free="<<free<<endl;
7104 if(displacement>0) mBaseRotationAmplitude*=0.1/displacement;
7105 if(mBaseRotationAmplitude<(0.02*M_PI/20.))
7107 mBaseRotationAmplitude=0.02*M_PI/20.;
7111 if(mBaseRotationAmplitude>(0.02*M_PI*20.))
7113 mBaseRotationAmplitude=0.02*M_PI*20.;
7121 VFN_DEBUG_EXIT(
"Molecule::TuneGlobalOptimRotationAmplitude()",5)
7127 if( (mClockFlipGroup>mClockConnectivityTable)
7128 &&(mClockFlipGroup>mClockRigidGroup))
return;
7129 VFN_DEBUG_ENTRY(
"Molecule::BuildFlipGroup()",5)
7130 TAU_PROFILE(
"Molecule::BuildFlipGroup()",
"void ()",TAU_DEFAULT);
7133 for(vector<MolAtom*>::const_iterator atom0=this->GetAtomList().begin();
7134 atom0!=this->GetAtomList().end();++atom0)
7137 if(pConn->size()<3)
continue;
7139 for(set<MolAtom*>::const_iterator pos1=pConn->begin();pos1!=pConn->end();++pos1)
7141 for(set<MolAtom*>::const_iterator pos2=pos1;pos2!=pConn->end();++pos2)
7143 if(*pos2==*pos1)
continue;
7147 bool foundRing=
false;
7148 for(set<MolAtom*>::const_iterator pos=pConn->begin();pos!=pConn->end();++pos)
7150 if((pos==pos1)||(pos==pos2))
continue;
7152 make_pair(*pos,set<MolAtom*>()));
7153 mvFlipGroup.back().mvRotatedChainList.back().second.insert(*atom0);
7155 mvFlipGroup.back().mvRotatedChainList.back().second);
7156 mvFlipGroup.back().mvRotatedChainList.back().second.erase(*atom0);
7157 set<MolAtom*>::const_iterator ringdetect1,ringdetect2;
7158 ringdetect1=find(
mvFlipGroup.back().mvRotatedChainList.back().second.begin(),
7159 mvFlipGroup.back().mvRotatedChainList.back().second.end(),
7161 ringdetect2=find(
mvFlipGroup.back().mvRotatedChainList.back().second.begin(),
7162 mvFlipGroup.back().mvRotatedChainList.back().second.end(),
7164 if( (ringdetect1!=
mvFlipGroup.back().mvRotatedChainList.back().second.end())
7165 ||(ringdetect2!=
mvFlipGroup.back().mvRotatedChainList.back().second.end()))
7168 unsigned long flipSize=0;
7169 for(list<pair<
const MolAtom *,set<MolAtom*> > >::const_iterator
7170 chain=
mvFlipGroup.back().mvRotatedChainList.begin();
7171 chain!=
mvFlipGroup.back().mvRotatedChainList.end();++chain)
7172 flipSize+=chain->second.size();
7179 make_pair(*atom0,set<MolAtom*>()));
7180 mvFlipGroup.back().mvRotatedChainList.back().second.insert(*atom0);
7182 mvFlipGroup.back().mvRotatedChainList.back().second);
7184 mvFlipGroup.back().mvRotatedChainList.back().second);
7185 mvFlipGroup.back().mvRotatedChainList.back().second.erase(*atom0);
7195 set<MolAtom*> fullset;
7196 for(list<pair<
const MolAtom *,set<MolAtom *> > >::iterator chain=pos->mvRotatedChainList.begin();
7197 chain!=pos->mvRotatedChainList.end();++chain)
7198 for(set<MolAtom *>::const_iterator at=chain->second.begin();at!=chain->second.end();++at)
7199 fullset.insert(*at);
7205 for(set<MolAtom *>::const_iterator at=fullset.begin();at!=fullset.end();++at)
7206 ct+=(*group)->count(*at);
7208 if((ct>0)&&(ct<(*group)->size())) {keep=
false;
break;}
7213 cout <<
"EXCLUDING flip group (breaking a rigid group): "
7214 <<pos->mpAtom0->GetName()<<
",exchanging bonds with "
7215 <<pos->mpAtom1->GetName()<<
" and "
7216 <<pos->mpAtom2->GetName()<<
", resulting in a 180deg rotation of atoms : ";
7217 for(set<MolAtom*>::iterator pos1=pos->mvRotatedChainList.begin()->second.begin();
7218 pos1!=pos->mvRotatedChainList.begin()->second.end();++pos1)
7219 cout<<(*pos1)->GetName()<<
" ";
7230 if(pos->mpAtom0->IsNonFlipAtom())
7232 cout <<
"EXCLUDING flip group (central atom is in the non-flip list): "
7233 <<pos->mpAtom0->GetName()<<
",exchanging bonds with "
7234 <<pos->mpAtom1->GetName()<<
" and "
7235 <<pos->mpAtom2->GetName()<<
", resulting in a 180deg rotation of atoms : ";
7236 for(set<MolAtom*>::iterator pos1=pos->mvRotatedChainList.begin()->second.begin();
7237 pos1!=pos->mvRotatedChainList.begin()->second.end();++pos1)
7238 cout<<(*pos1)->GetName()<<
" ";
7252 for(list<FlipGroup>::iterator pos=
mvFlipGroup.begin();
7255 if(pos->mvRotatedChainList.begin()->first==pos->mpAtom0)
7257 cout <<
"Flip group from atom "
7258 <<pos->mpAtom0->GetName()<<
",exchanging bonds with "
7259 <<pos->mpAtom1->GetName()<<
" and "
7260 <<pos->mpAtom2->GetName()<<
", resulting in a 180� rotation of atoms : ";
7261 for(set<MolAtom*>::iterator pos1=pos->mvRotatedChainList.begin()->second.begin();
7262 pos1!=pos->mvRotatedChainList.begin()->second.end();++pos1)
7263 cout<<(*pos1)->GetName()<<
" ";
7267 cout <<
"Flip group with respect to: "
7268 <<pos->mpAtom1->GetName()<<
"-"
7269 <<pos->mpAtom0->GetName()<<
"-"
7270 <<pos->mpAtom2->GetName()<<
" : ";
7271 for(list<pair<
const MolAtom *,set<MolAtom*> > >::const_iterator
7272 chain=pos->mvRotatedChainList.begin();
7273 chain!=pos->mvRotatedChainList.end();++chain)
7275 cout<<
" -"<<chain->first->GetName()<<
":";
7276 for(set<MolAtom*>::const_iterator pos1=chain->second.begin();
7277 pos1!=chain->second.end();++pos1)
7278 cout<<(*pos1)->GetName()<<
" ";
7290 cout <<
" -> NOT a free flip, d(llk)="<<dllk;
7293 else cout <<
" -> free flip, d(llk)="<<dllk;
7299 mClockFlipGroup.
Click();
7300 VFN_DEBUG_EXIT(
"Molecule::BuildFlipGroup()",5)
7306 map<const MolBond*, set<const StretchMode*> > vpBond;
7310 if(mode->mpBond!=0) vpBond[mode->mpBond].insert(&(*mode));
7311 for(map<const MolBond*,REAL>::const_iterator pos=mode->mvpBrokenBond.begin();
7312 pos!=mode->mvpBrokenBond.end();++pos)
7313 vpBond[pos->first].insert(&(*mode));
7317 for(map<const MolBond*,REAL>::const_iterator pos=mode->mvpBrokenBond.begin();
7318 pos!=mode->mvpBrokenBond.end();++pos)
7319 vpBond[pos->first].insert(&(*mode));
7323 for(map<const MolBond*,REAL>::const_iterator pos=mode->mvpBrokenBond.begin();
7324 pos!=mode->mvpBrokenBond.end();++pos)
7325 vpBond[pos->first].insert(&(*mode));
7328 for(map<const MolBond*,REAL>::const_iterator pos=mode->mvpBrokenBond.begin();
7329 pos!=mode->mvpBrokenBond.end();++pos)
7330 vpBond[pos->first].insert(&(*mode));
7334 map<const MolBondAngle*, set<const StretchMode*> > vpAngle;
7337 for(map<const MolBondAngle*,REAL>::const_iterator pos=mode->mvpBrokenBondAngle.begin();
7338 pos!=mode->mvpBrokenBondAngle.end();++pos)
7339 vpAngle[pos->first].insert(&(*mode));
7344 if(mode->mpBondAngle!=0) vpAngle[mode->mpBondAngle].insert(&(*mode));
7345 for(map<const MolBondAngle*,REAL>::const_iterator pos=mode->mvpBrokenBondAngle.begin();
7346 pos!=mode->mvpBrokenBondAngle.end();++pos)
7347 vpAngle[pos->first].insert(&(*mode));
7351 for(map<const MolBondAngle*,REAL>::const_iterator pos=mode->mvpBrokenBondAngle.begin();
7352 pos!=mode->mvpBrokenBondAngle.end();++pos)
7353 vpAngle[pos->first].insert(&(*mode));
7356 for(map<const MolBondAngle*,REAL>::const_iterator pos=mode->mvpBrokenBondAngle.begin();
7357 pos!=mode->mvpBrokenBondAngle.end();++pos)
7358 vpAngle[pos->first].insert(&(*mode));
7362 map<const MolDihedralAngle*,set<const StretchMode*> > vpDihed;
7365 for(map<const MolDihedralAngle*,REAL>::const_iterator pos=mode->mvpBrokenDihedralAngle.begin();
7366 pos!=mode->mvpBrokenDihedralAngle.end();++pos)
7367 vpDihed[pos->first].insert(&(*mode));
7371 for(map<const MolDihedralAngle*,REAL>::const_iterator pos=mode->mvpBrokenDihedralAngle.begin();
7372 pos!=mode->mvpBrokenDihedralAngle.end();++pos)
7373 vpDihed[pos->first].insert(&(*mode));
7378 if(mode->mpDihedralAngle!=0) vpDihed[mode->mpDihedralAngle].insert(&(*mode));
7379 for(map<const MolDihedralAngle*,REAL>::const_iterator pos=mode->mvpBrokenDihedralAngle.begin();
7380 pos!=mode->mvpBrokenDihedralAngle.end();++pos)
7381 vpDihed[pos->first].insert(&(*mode));
7385 for(map<const MolDihedralAngle*,REAL>::const_iterator pos=mode->mvpBrokenDihedralAngle.begin();
7386 pos!=mode->mvpBrokenDihedralAngle.end();++pos)
7387 vpDihed[pos->first].insert(&(*mode));
7389 for(map<
const MolBond*,set<const StretchMode*> >::const_iterator pos=vpBond.begin();pos!=vpBond.end();++pos)
7391 cout<<
"Bond "<<pos->first->GetName()<<
" is modified by the stretch modes:"<<endl;
7392 for(set<const StretchMode*>::const_iterator mode=pos->second.begin();mode!=pos->second.end();++mode)
7395 (*mode)->Print(cout);
7399 for(map<
const MolBondAngle*,set<const StretchMode*> >::const_iterator pos=vpAngle.begin();pos!=vpAngle.end();++pos)
7401 cout<<
"Bond Angle "<<pos->first->GetName()<<
" is modified by the stretch modes:"<<endl;
7402 for(set<const StretchMode*>::const_iterator mode=pos->second.begin();mode!=pos->second.end();++mode)
7405 (*mode)->Print(cout);
7409 for(map<
const MolDihedralAngle*,set<const StretchMode*> >::const_iterator pos=vpDihed.begin();pos!=vpDihed.end();++pos)
7411 cout<<
"Dihedral Angle "<<pos->first->GetName()<<
" is modified by the stretch modes:"<<endl;
7412 for(set<const StretchMode*>::const_iterator mode=pos->second.begin();mode!=pos->second.end();++mode)
7415 (*mode)->Print(cout);
7426 int nb=mode->mvpBrokenDihedralAngle.size()+mode->mvpBrokenBondAngle.size()+mode->mvpBrokenBond.size();
7427 if(mode->mpBond!=0) nb -= 1;
7435 int nb=mode->mvpBrokenDihedralAngle.size()+mode->mvpBrokenBondAngle.size()+mode->mvpBrokenBond.size();
7436 if(mode->mpBondAngle!=0) nb -= 1;
7444 int nb=mode->mvpBrokenDihedralAngle.size()+mode->mvpBrokenBondAngle.size()+mode->mvpBrokenBond.size();
7445 if(mode->mpDihedralAngle!=0) nb -= 1;
7452 if(mode->mvpBrokenDihedralAngle.size()+mode->mvpBrokenBondAngle.size()+mode->mvpBrokenBond.size()==0)
7462 map<MolAtom*,set<MolAtom*> > vBoundAtoms;
7464 for(vector<MolAtom*>::iterator pos=this->GetAtomList().begin();pos!=this->GetAtomList().end();++pos)
7467 for(vector<MolAtom*>::iterator pat1=this->GetAtomList().begin();pat1!=this->GetAtomList().end();++pat1)
7469 vBoundAtoms[*pat1]=set0;
7470 for(vector<MolAtom*>::iterator pat2=this->GetAtomList().begin();pat2!=this->GetAtomList().end();++pat2)
7474 for(list<StretchModeBondLength>::iterator pstretch=this->GetStretchModeBondLengthList().begin();
7475 pstretch!=this->GetStretchModeBondLengthList().end();++pstretch)
7477 set<MolAtom *>::iterator pos1=pstretch->mvTranslatedAtomList.find(*pat1),
7478 pos2=pstretch->mvTranslatedAtomList.find(*pat2);
7479 if( ((pos1==pstretch->mvTranslatedAtomList.end())&&(pos2!=pstretch->mvTranslatedAtomList.end()))
7480 ||((pos1!=pstretch->mvTranslatedAtomList.end())&&(pos2==pstretch->mvTranslatedAtomList.end())))
7482 vBoundAtoms[*pat1].erase(*pat2);
7491 for(list<StretchModeBondAngle>::iterator pstretch=this->GetStretchModeBondAngleList().begin();
7492 pstretch!=this->GetStretchModeBondAngleList().end();++pstretch)
7495 if((*pat1==pstretch->mpAtom1)||(*pat2==pstretch->mpAtom1))
continue;
7497 set<MolAtom *>::iterator pos1=pstretch->mvRotatedAtomList.find(*pat1),
7498 pos2=pstretch->mvRotatedAtomList.find(*pat2);
7500 if( ((pos1==pstretch->mvRotatedAtomList.end())&&(pos2!=pstretch->mvRotatedAtomList.end()))
7501 ||((pos1!=pstretch->mvRotatedAtomList.end())&&(pos2==pstretch->mvRotatedAtomList.end())))
7503 vBoundAtoms[*pat1].erase(*pat2);
7512 for(list<StretchModeTorsion>::iterator pstretch=this->GetStretchModeTorsionList().begin();
7513 pstretch!=this->GetStretchModeTorsionList().end();++pstretch)
7515 set<MolAtom *>::iterator pos1=pstretch->mvRotatedAtomList.find(*pat1),
7516 pos2=pstretch->mvRotatedAtomList.find(*pat2);
7518 if( (pos1!=pstretch->mvRotatedAtomList.end())&&(pos2==pstretch->mvRotatedAtomList.end())
7519 &&(*pat2!=pstretch->mpAtom1) &&(*pat2!=pstretch->mpAtom2) )
7521 vBoundAtoms[*pat1].erase(*pat2);
7526 if( (pos1==pstretch->mvRotatedAtomList.end())&&(pos2!=pstretch->mvRotatedAtomList.end())
7527 &&(*pat1!=pstretch->mpAtom1) && (*pat1!=pstretch->mpAtom2) )
7529 vBoundAtoms[*pat1].erase(*pat2);
7539 set<set<MolAtom*> > vBoundGroups;
7540 for(map<
MolAtom*,set<MolAtom*> >::iterator pos=vBoundAtoms.begin();pos!=vBoundAtoms.end();++pos)
7543 cout<<
"Non-flexible group from "<<pos->first->GetName()<<
": ";
7544 for(set<MolAtom*>::const_iterator atom=pos->second.begin();atom!=pos->second.end();++atom)
7545 cout<<(*atom)->GetName()<<
",";
7550 for(set<MolAtom *>::iterator at=(*pr)->begin();at!=(*pr)->end();++at)
7551 pos->second.erase(*at);
7552 if(pos->second.size()>1) vBoundGroups.insert(pos->second);
7555 for(set<set<MolAtom*> >::iterator pos=vBoundGroups.begin();pos!=vBoundGroups.end();++pos)
7557 cout<<
"Non-flexible group:";
7558 for(set<MolAtom*>::const_iterator atom=pos->begin();atom!=pos->end();++atom)
7559 cout<<(*atom)->GetName()<<
",";
7565 for(set<set<MolAtom*> >::iterator pos=vBoundGroups.begin();pos!=vBoundGroups.end();++pos)
7568 for(vector<MolBond*>::iterator pr=this->GetBondList().begin();pr!=this->GetBondList().end();++pr)
7569 if( (pos->find(&(*pr)->GetAtom1())!=pos->end())
7570 ||(pos->find(&(*pr)->GetAtom2())!=pos->end())) vb.insert(*pr);
7572 set<MolBondAngle*> va;
7573 for(vector<MolBondAngle*>::iterator pr=this->GetBondAngleList().begin();pr!=this->GetBondAngleList().end();++pr)
7574 if( (pos->find(&(*pr)->GetAtom1())!=pos->end())
7575 ||(pos->find(&(*pr)->GetAtom2())!=pos->end())
7576 ||(pos->find(&(*pr)->GetAtom3())!=pos->end())) va.insert(*pr);
7578 set<MolDihedralAngle*> vd;
7579 for(vector<MolDihedralAngle*>::iterator pr=this->GetDihedralAngleList().begin();pr!=this->GetDihedralAngleList().end();++pr)
7580 if( (pos->find(&(*pr)->GetAtom1())!=pos->end())
7581 ||(pos->find(&(*pr)->GetAtom2())!=pos->end())
7582 ||(pos->find(&(*pr)->GetAtom3())!=pos->end())
7583 ||(pos->find(&(*pr)->GetAtom4())!=pos->end())) vd.insert(*pr);
7585 set<MolAtom*> tmp= *pos;
7594 for(vector<MolAtom*>::iterator at=this->GetAtomList().begin();at!=this->GetAtomList().end();++at)
7597 for(set<MolAtom *>::iterator at=(*pr)->begin();at!=(*pr)->end();++at)
7599 cout<<
"Full MD atom group:"<<endl<<
" ";
7601 cout<<(*pos)->GetName()<<
" ";
7607 for(set<MolAtom*>::const_iterator at=pos->mvpAtom.begin();at!=pos->mvpAtom.end();++at)
7609 cout<<
"Full MD atom group:"<<endl<<
" ";
7611 cout<<(*pos)->GetName()<<
" ";
7614 mClockMDAtomGroup.
Click();
7624 VFN_DEBUG_ENTRY(
"Molecule::UpdateScattCompList()",5)
7625 TAU_PROFILE(
"Molecule::UpdateScattCompList()",
"void ()",TAU_DEFAULT);
7628 for(
long i=0;i<nb;++i)
7638 #ifdef RIGID_BODY_STRICT_EXPERIMENTAL
7644 (*pos)->mQuat.Normalize();
7646 REAL x0=0,y0=0,z0=0;
7647 for(set<unsigned int>::iterator at=(*pos)->mvIdx.begin();at!=(*pos)->mvIdx.end();++at)
7658 for(set<unsigned int>::iterator at=(*pos)->mvIdx.begin();at!=(*pos)->mvIdx.end();++at)
7661 (*pos)->mQuat.RotateVector(x,y,z);
7670 REAL x0=0,y0=0,z0=0;
7673 for(
long i=0;i<nb;++i)
7689 for(
long i=0;i<nb;++i)
7698 for(
long i=0;i<nb;++i)
7704 for(
long i=0;i<nb;++i)
7711 for(
long i=0;i<nb;++i)
7718 VFN_DEBUG_EXIT(
"Molecule::UpdateScattCompList()",5)
7722 VFN_DEBUG_ENTRY(
"Molecule::FindAtom():"<<name,4)
7723 vector<MolAtom*>::reverse_iterator rpos;
7725 if(name==(*rpos)->GetName())
7727 VFN_DEBUG_EXIT(
"Molecule::FindAtom():"<<name<<
"...FOUND !",4)
7730 VFN_DEBUG_EXIT(
"Molecule::FindAtom():"<<name<<
"...NOT FOUND !",4)
7735 vector<MolAtom*>::const_reverse_iterator rpos;
7738 if(name==(*rpos)->GetName())
return rpos;
7743 VFN_DEBUG_ENTRY(
"Molecule::InitOptions",7)
7744 static string Flexname;
7745 static string Flexchoices[3];
7747 static string FlipName;
7748 static string FlipChoice[2];
7750 static string autoOptimizeConformationName;
7751 static string autoOptimizeConformationChoices[2];
7753 static string optimizeOrientationName;
7754 static string optimizeOrientationChoices[2];
7756 static string moleculeCenterName;
7757 static string moleculeCenterChoices[2];
7759 static bool needInitNames=
true;
7760 if(
true==needInitNames)
7762 Flexname=
"Flexibility Model";
7763 Flexchoices[0]=
"Automatic from Restraints, relaxed - RECOMMENDED";
7764 Flexchoices[1]=
"Rigid Body";
7765 Flexchoices[2]=
"Automatic from Restraints, strict";
7768 FlipName=
"Enable Flipping";
7769 FlipChoice[0]=
"Yes";
7772 autoOptimizeConformationName=
"Auto Optimize Starting Conformation";
7773 autoOptimizeConformationChoices[0]=
"Yes";
7774 autoOptimizeConformationChoices[1]=
"No";
7776 optimizeOrientationName=
"Optimize Orientation";
7777 optimizeOrientationChoices[0]=
"Yes";
7778 optimizeOrientationChoices[1]=
"No";
7780 moleculeCenterName=
"Rotation Center";
7781 moleculeCenterChoices[0]=
"Geometrical center (recommended)";
7782 moleculeCenterChoices[1]=
"User-chosen Atom";
7784 needInitNames=
false;
7790 mFlipModel.Init(2, &FlipName, FlipChoice);
7791 mFlipModel.SetChoice(0);
7795 autoOptimizeConformationChoices);
7804 VFN_DEBUG_EXIT(
"Molecule::InitOptions",7)
7808 mpAtom0(&at0),mpAtom1(&at1),mpAtom2(&at2),mNbTest(0),mNbAccept(0)
7814 TAU_PROFILE(
"Molecule::FlipAtomGroup(FlipGroup&)",
"void (...)",TAU_DEFAULT);
7828 const REAL norm01=sqrt(v01x*v01x+v01y*v01y+v01z*v01z+1e-7);
7829 v01x /= norm01;v01y /= norm01;v01z /= norm01;
7834 const REAL norm02=sqrt(v02x*v02x+v02y*v02y+v02z*v02z+1e-7);
7835 v02x /= norm02;v02y /= norm02;v02z /= norm02;
7840 const REAL norm12=sqrt(v12x*v12x+v12y*v12y+v12z*v12z+1e-7);
7841 v12x /= norm12;v12y /= norm12;v12z /= norm12;
7846 const REAL norm0m=sqrt(v0mx*v0mx+v0my*v0my+v0mz*v0mz+1e-7);
7847 v0mx /= norm0m;v0my /= norm0m;v0mz /= norm0m;
7849 if(fabs(v01x*v02x+v01y*v02y+v01z*v02z)
7850 >0.05*sqrt(abs( (v01x*v01x+v01y*v01y+v01z*v01z)
7851 *(v02x*v02x+v02y*v02y+v02z*v02z))))
7853 REAL v012x=v01y*v02z-v01z*v02y;
7854 REAL v012y=v01z*v02x-v01x*v02z;
7855 REAL v012z=v01x*v02y-v01y*v02x;
7856 const REAL norm012=sqrt(v012x*v012x+v012y*v012y+v012z*v012z+1e-7);
7857 v012x /= norm012;v012y /= norm012;v012z /= norm012;
7860 for(list<pair<
const MolAtom *,set<MolAtom*> > >::const_iterator
7864 REAL v03x=chain->first->X()-group.
mpAtom0->X();
7865 REAL v03y=chain->first->Y()-group.
mpAtom0->Y();
7866 REAL v03z=chain->first->Z()-group.
mpAtom0->Z();
7867 const REAL norm03=sqrt( v03x*v03x + v03y*v03y + v03z*v03z +1e-7);
7868 v03x /= norm03;v03y /= norm03;v03z /= norm03;
7870 const REAL a1=v012x*v03x+v012y*v03y+v012z*v03z;
7871 const REAL a2= v0mx*v03x+ v0my*v03y+ v0mz*v03z;
7872 const REAL a3= v12x*v03x+ v12y*v03y+ v12z*v03z;
7873 REAL angle = -a1/sqrt(1-a3*a3+1e-7);
7882 else angle = asin(angle);
7884 if(a2<0) angle=M_PI-angle;
7886 chain->second,2*angle,keepCenter);
7898 #ifdef RIGID_BODY_STRICT_EXPERIMENTAL
7903 (*pos)->mQuat.Normalize();
7905 REAL x0=0,y0=0,z0=0;
7906 for(set<MolAtom *>::iterator at=(*pos)->begin();at!=(*pos)->end();++at)
7917 for(set<MolAtom *>::iterator at=(*pos)->begin();at!=(*pos)->end();++at)
7919 REAL x=(*at)->GetX()-x0, y=(*at)->GetY()-y0, z=(*at)->GetZ()-z0;
7920 (*pos)->mQuat.RotateVector(x,y,z);
7921 (*at)->SetX(x+x0+(*pos)->mX);
7922 (*at)->SetY(y+y0+(*pos)->mY);
7923 (*at)->SetZ(z+z0+(*pos)->mZ);
7930 (*pos)->mQuat.Q0()=1;
7931 (*pos)->mQuat.Q1()=0;
7932 (*pos)->mQuat.Q2()=0;
7933 (*pos)->mQuat.Q3()=0;
7938 #ifdef __WX__CRYST__
7941 VFN_DEBUG_ENTRY(
"Molecule::WXCreate()",5)
7943 VFN_DEBUG_EXIT("
Molecule::WXCreate()",5)
7944 return mpWXCrystObj;
7951 VFN_DEBUG_ENTRY(
"ZScatterer2Molecule()",6)
7954 REAL x0=0,y0=0,z0=0;
7955 for(
unsigned int i=0;i<nb;++i)
7970 .GetObj(i).GetZBondLength())));
7971 if(pLength->IsFixed())
7973 pLength->
GetValue(),.01,.02,
false);
7975 if(pLength->IsLimited())
7977 (pLength->
GetMin()+pLength->
GetMax())/2.,.01,.02,
false);
7982 .GetObj(i).GetZAngle())));
7983 if(pAngle->IsFixed())
7988 if(pAngle->IsLimited())
7996 .GetObj(i).GetZDihedralAngle())));
7997 MolAtom *p1=&(mol->GetAtom(i));
8006 if(pDihed->IsFixed())
8009 if(((pDihed->
GetMax()-pDihed->
GetMax())<0.3)&&(i>2)&&(pDihed->IsLimited()))
8015 mol->GetAtom(i).SetOccupancy(scatt->
GetZAtomRegistry().GetObj(i).GetOccupancy());
8018 CrystVector_REAL x(nb),y(nb),z(nb),radius(nb);
8019 vector<pair<const ScatteringPowerAtom *,long> > scattpow(nb);
8020 for(
unsigned int i=0;i<nb;++i)
8022 x(i)=mol->GetAtom(i).GetX();
8023 y(i)=mol->GetAtom(i).GetY();
8024 z(i)=mol->GetAtom(i).GetZ();
8028 scattpow[i].first=0;
8032 radius(i)=mol->GetAtom(i).GetScatteringPower().
GetRadius();
8034 (&(mol->GetAtom(i).GetScatteringPower()));
8038 for(
unsigned int i=0;i<nb;++i)
8040 if(scattpow[i].first==0)
continue;
8041 const REAL x1=x(i),y1=y(i),z1=z(i);
8045 for(
unsigned int j=i+1;j<nb;++j)
8047 if(scattpow[j].first==0)
continue;
8048 const REAL dist=sqrt(x(j)*x(j)+y(j)*y(j)+z(j)*z(j));
8050 if(dist<(1.10*(radius(i)+radius(j))))
8052 if((1!=scattpow[i].second)||(1!=scattpow[j].second))
8054 mol->
AddBond(mol->GetAtom(i),mol->GetAtom(j),dist,.01,.02,
false);
8066 for(set<MolAtom*>::const_iterator pos1=pos->second.begin();
8067 pos1!=pos->second.end();++pos1)
8069 for(set<MolAtom*>::const_iterator pos2=pos1;
8070 pos2!=pos->second.end();++pos2)
8072 if(pos2==pos1)
continue;
8073 if(mol->
FindBondAngle(**pos1,*(pos->first),**pos2)== mol->GetBondAngleList().end())
8075 GetBondAngle(**pos1,*(pos->first),**pos2),0.01,0.02,
false);
8087 VFN_DEBUG_EXIT(
"ZScatterer2Molecule()",6)
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++.
REAL GetBondLength(const MolAtom &at1, const MolAtom &at2)
Get The Bond Length between two atoms.
float string2floatC(const string &s)
Function to convert a substring to a floating point value, imposing a C locale (using '.
Molecule * ZScatterer2Molecule(ZScatterer *scatt)
Converter from ZScatterer to a Molecule object.
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.
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 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.
REAL GetDihedralAngle(const MolAtom &at1, const MolAtom &at2, const MolAtom &at3, const MolAtom &at4)
Get The dihedral angle defined by 4 atoms.
REAL GetBondAngle(const MolAtom &at1, const MolAtom &at2, const MolAtom &at3)
Get The Bond Angle of 3 atoms.
ObjRegistry< Crystal > gCrystalRegistry("List of all Crystals")
Global registry for all Crystal objects.
Crystal class: Unit cell, spacegroup, scatterers.
ScatteringPower & GetScatteringPower(const string &name)
Find a ScatteringPower from its name. Names must be unique in a given Crystal.
REAL GetDynPopCorr(const Scatterer *pscatt, unsigned int component) const
Access the Dynamical Occupancy Correction for a given component (atom) in a given Scatterer.
Exception class for ObjCryst++ library.
Class to store POV-Ray output options.
bool mShowHydrogens
Show hydrogens ?
REAL mXmin
Display limits in reduced coordinates.
bool mShowLabel
Show labels ?
Structure holding 3 coordinates, or deriviatives with respect to each of these coordinates.
MolAtom : atom inside a Molecule.
virtual ~MolAtom()
Destructor.
const ScatteringPower * mpScattPow
ScatteringPower.
REAL mX
Cartesian oordinates in the Molecule reference frame.
size_t int_ptr() const
Access to the integer address of this object, for unique identification from python.
bool mIsInRing
Is the atom in a ring ?
Molecule * mpMol
Parent Molecule.
bool IsDummy() const
Returns true if this is a dummy atom, i.e.
bool mIsNonFlipAtom
Can the atom be flipped (this is used for optically active atom which should keep their absolute conf...
void SetX(const REAL) const
Set the X,Y,Z coordinate - this is const because sometimes their coordinate must be changed even thou...
string mName
Name for this atom.
bool IsNonFlipAtom() const
Can this atom be flipped (return=false) or should its absolute configuration be kept (return=true)
REAL mOccupancy
Occupancy.
void SetNonFlipAtom(const bool nonflip)
Set a flag to prevent the atom's absolute configuration to be changed.
void SetIsInRing(const bool r) const
Flag this atom as being in a ring (or not).
Bond between two atoms, also a restraint on the associated bond length.
size_t int_ptr() const
Access to the integer address of this object, for unique identification from python.
virtual ~MolBond()
Destructor.
REAL mLLK
Stored log(likelihood)
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...
virtual REAL GetLogLikelihood() const
Get -ln(likelihood) for this restraint.
string GetName() const
Name of the bond, e.g. "C3-O4".
MolBond(MolAtom &atom1, MolAtom &atom2, const REAL length, const REAL sigma, const REAL delta, Molecule &parent, const REAL bondOrder=1.)
Constructor.
void CalcGradient(std::map< MolAtom *, XYZ > &m) const
Calc log(likelihood) gradient - versus all atomic coordinates.
Molecule * mpMol
Parent Molecule.
REAL mDerivLLKCoeff
The factor used to change the derivative of the length/angle, to the derivative of the log(likelihood...
XYZ mDerivAtom1
Derivatives of the bond length with respect to the coordinates of the atoms.
Bond angle restraint between 3 atoms.
virtual REAL GetLogLikelihood() const
Get -ln(likelihood) for this restraint.
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...
XYZ mDerivAtom1
Partial derivatives of the angle with respect to the coordinates of the atoms.
void CalcGradient(std::map< MolAtom *, XYZ > &m) const
Calc log(likelihood) gradient - versus all atomic coordinates.
MolBondAngle(MolAtom &atom1, MolAtom &atom2, MolAtom &atom3, const REAL angle, const REAL sigma, const REAL delta, Molecule &parent)
Constructor.
size_t int_ptr() const
Access to the integer address of this object, for unique identification from python.
REAL mLLK
Stored log(likelihood)
Molecule * mpMol
Parent Molecule.
vector< MolAtom * > mvpAtom
The vector of the 3 atoms involved in the bond angle.
REAL mDerivLLKCoeff
The factor used to change the derivative of the length/angle, to the derivative of the log(likelihood...
virtual ~MolBondAngle()
Destructor.
Dihedral angle restraint between 4 atoms.
XYZ mDerivAtom1
Partial derivatives of the angle with respect to the coordinates of the atoms.
size_t int_ptr() const
Access to the integer address of this object, for unique identification from python.
vector< MolAtom * > mvpAtom
The vector of the 4 atoms involved in the bond angle.
REAL mDerivLLKCoeff
The factor used to change the derivative of the length/angle, to the derivative of the log(likelihood...
virtual ~MolDihedralAngle()
Destructor.
void CalcGradient(std::map< MolAtom *, XYZ > &m) const
Calc log(likelihood) gradient - versus all atomic coordinates.
REAL mLLK
Stored log(likelihood)
MolDihedralAngle(MolAtom &atom1, MolAtom &atom2, MolAtom &atom3, MolAtom &atom4, const REAL angle, const REAL sigma, const REAL delta, Molecule &parent)
Constructor.
Molecule * mpMol
Parent Molecule.
virtual REAL GetLogLikelihood() const
Get -ln(likelihood) for this restraint.
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...
size_t int_ptr() const
Access to the integer address of this object, for unique identification from python.
A quaternion class, used to represent the orientation of the molecule.
Quaternion GetConjugate() const
Get the conjugate of this quaternion (== the inverse if unit quaternion)
Quaternion operator*(const Quaternion &q) const
Quaternion multiplication.
REAL mQ0
The components of the quaternion z=(q0,v) with v=(q1,q2,q3)
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.
void Normalize() const
Re-normalize the quaternion to unity.
Quaternion()
Default constructor, yields q=(1,0,0,0)
void RotateVector(REAL &v1, REAL &v2, REAL &v3) const
Rotate vector v=(v1,v2,v3). The rotated components are directly written.
Rigid groups of atoms inside a molecule.
Quaternion mQuat
The unit quaternion defining the orientation - this is used during optimizations to rotate all atoms ...
size_t int_ptr() const
Access to the integer address of this object, for unique identification from python.
REAL mX
The translation of all the atoms as a group The values will be resetted whenever entering or leaving ...
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...
Molecule * mpMol
The Molecule corresponding to this stretch mode.
std::map< const MolAtom *, XYZ > mDerivXYZ
Derivative of the atomic positions versus a change of the bond length.
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...
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...
REAL mBaseAmplitude
The recommended change amplitude, for a base global optimization displacement, to obtain an average 0...
REAL mLLKDeriv
Derivative of the Molecule's Log(likelihood) versus a change of the bond length.
Group of atoms for random moves changing a bond length.
virtual void RandomStretch(const REAL amplitude, const bool keepCenter=true)
Move the atoms according to this mode, randomly.
MolAtom * mpAtom0
The first atom (fixed).
const MolBond * mpBond
The (optional) bond length which this stretch mode should respect.
StretchModeBondLength(MolAtom &at0, MolAtom &at1, const MolBond *pBond)
Constructor If pBond!=0, the bond length restraint is respected.
virtual void Print(ostream &os, bool full=true) const
Print one-line list of atoms moved.
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...
set< MolAtom * > mvTranslatedAtomList
The set of atoms that are to be translated, including at1.
MolAtom * mpAtom1
The second atom (first atom moved)
virtual void Stretch(const REAL change, const bool keepCenter=true)
Move the atoms according to this mode.
Atoms moved when changing a bond angle.
MolAtom * mpAtom2
The third atom.
set< MolAtom * > mvRotatedAtomList
The set of atoms that are to be rotated around the direction going through at1 and perpendicular to t...
virtual void Print(ostream &os, bool full=true) const
Print one-line list of atoms moved.
virtual void RandomStretch(const REAL amplitude, const bool keepCenter=true)
Move the atoms according to this mode, randomly.
MolAtom * mpAtom0
The first atom.
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 void Stretch(const REAL change, const bool keepCenter=true)
Move the atoms according to this mode.
StretchModeBondAngle(MolAtom &at0, MolAtom &at1, MolAtom &at2, const MolBondAngle *pBondAngle)
Constructor If pBondAngle!=0, the bond angle length restraint is respected.
MolAtom * mpAtom1
The second atom.
const MolBondAngle * mpBondAngle
The (optional) bond angle restraint which this stretch mode should respect.
Atoms moved when rotated around a bond at0-at1-at2-at3.
MolAtom * mpAtom1
The first atom.
MolAtom * mpAtom2
The second atom.
StretchModeTorsion(MolAtom &at1, MolAtom &at2, const MolDihedralAngle *pDihedralAngle)
Constructor If pDihedralAngle!=0, the dihedral angle length restraint is respected.
virtual void Print(ostream &os, bool full=true) const
Print one-line list of atoms moved.
const MolDihedralAngle * mpDihedralAngle
The (optional) bond angle restraint which this stretch mode should respect.
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.
virtual void Stretch(const REAL change, const bool keepCenter=true)
Move the atoms according to this mode.
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...
Atoms moved between two other atoms, using a "twist" of their positions - only small twists of their ...
MolAtom * mpAtom2
The second atom.
virtual void Print(ostream &os, bool full=true) const
Print one-line list of atoms moved.
set< MolAtom * > mvRotatedAtomList
The set of atoms that are to be rotated around at1-at2.
StretchModeTwist(MolAtom &at1, MolAtom &at2)
Constructor If pDihedralAngle!=0, the dihedral angle length restraint is respected.
virtual void Stretch(const REAL change, const bool keepCenter=true)
Move the atoms according to this mode.
MolAtom * mpAtom1
The first atom.
virtual void RandomStretch(const REAL amplitude, const bool keepCenter=true)
Move the atoms according to this mode, randomly.
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...
Groups of atoms that can be moved using molecular dynamics principles, taking a list of restraints as...
MDAtomGroup()
Default constructor.
void Print(ostream &os, bool full=true) const
Print one-line list of atoms moved.
Molecule : class for complex scatterer descriptions using cartesian coordinates with bond length/angl...
virtual const ScatteringComponentList & GetScatteringComponentList() const
Get the list of all scattering components for this scatterer.
vector< MolBondAngle * >::iterator RemoveBondAngle(const MolBondAngle &, const bool del=true)
Remove a BondAngle.
void BuildMDAtomGroups()
Find groups of atoms that cannot be moved relatively to each other using the free or non-free stretch...
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.
map< MolAtom *, set< MolAtom * > > mConnectivityTable
Connectivity table: for each atom, keep the list of atoms bonded to it.
void InitOptions()
Build options for this object.
void BuildStretchModeTorsion()
Build the groups of atoms moved when changing a dihedral angle, while respecting the Molecule restrai...
void SetDeleteSubObjInDestructor(const bool b)
Set whether to delete the MolAtoms, MolBonds, MolBondAngles and MolDihedralAngles in the destructor.
vector< MolDihedralAngle * > mvpDihedralAngle
The list of dihedral angles.
std::list< StretchMode * > mvpStretchModeFree
Groups of StretchMode not breaking any restraint (unless the one they are associated to)
list< RotorGroup > mvRotorGroupTorsionSingleChain
List of RotorGroups corresponding to free torsion bonds, but with only one chain of atoms listed.
void SetCenterAtom(const MolAtom &at)
Get the atom defining the origin of the Molecule Equal to 0 if no atom as been set.
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.
const MolAtom * GetCenterAtom() const
Get the atom defining the origin of the Molecule Equal to 0 if no atom as been set.
virtual ostream & POVRayDescription(ostream &os, const CrystalPOVRayOptions &options) const
virtual void EndOptimization()
This should be called by any optimization class at the end of an optimization.
RefObjOpt mFlexModel
OPtion for the different types of flexibility possible for this molecule: rigid body,...
list< MolRing > mvRing
The list of rings.
virtual void SetName(const string &name)
Name of the object.
void BuildRotorGroup()
Build the groups of atoms that will be rotated during global optimization.
void BuildStretchModeBondAngle()
Build the groups of atoms moved when changing a bond angle, while respecting the Molecule restraints.
virtual unsigned int GetNbLSQFunction() const
Number of LSQ functions.
std::vector< RigidGroup * > mvRigidGroup
Rigid groups of atoms.
std::vector< MolZAtom > mAsZMatrix
The Molecule, as a lightweight ZMatrix, for export purposes.
void BuildConnectivityTable() const
Build the Connectivity table.
list< MDAtomGroup > mvMDAtomGroup
Groups of atoms that should be moved according to molecular dynamics principles.
RefObjOpt mMoleculeCenter
Option to choose the center of rotation of the Molecule for the global orientation either as the geom...
void FlipAtomGroup(const FlipGroup &, const bool keepCenter=true)
Flip a group of atom. See Molecule::FlipGroup.
list< StretchModeTwist > mvStretchModeTwist
List of StretchModeTwist.
void UpdateScattCompList() const
Update the Molecule::mScattCompList from the cartesian coordinates of all atoms, and the orientation ...
vector< MolDihedralAngle * >::iterator RemoveDihedralAngle(const MolDihedralAngle &, const bool del=true)
Remove a dihedral angle.
void RigidifyWithDihedralAngles()
Add dihedral angles so as to rigidify the Molecule.
list< StretchModeBondAngle > mvStretchModeBondAngle
List of StretchModeBondLength.
const std::vector< RigidGroup * > & GetRigidGroupList() const
List of rigid group of atoms.
virtual void InitRefParList()
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,...
bool mDeleteSubObjInDestructor
Base Rotation amplitude (in radians) for the Molecule, so that the average atomic displacement is equ...
REAL mMDMoveEnergy
Relative energy of molecule during molecular dynamics move Default: 40, 10 (slow conformation change)...
Quaternion mQuat
The unit quaternion defining the orientation.
CrystVector_REAL mLSQObs
Current LSQ Calc - one value for each restraint (bond distance, angle or dihedral angle ideal values)
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input From stream.
vector< MolAtom * >::iterator RemoveAtom(MolAtom &, const bool del=true)
Remove an atom.
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.
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,...
RefinableObjClock & GetRigidGroupClock()
Get the clock associated to the list of rigid groups (clicked also whenever a rigid group is modified...
virtual void TagNewBestConfig() const
During a global optimization, tells the object that the current config is the latest "best" config.
CrystVector_REAL mLSQWeight
Current LSQ Calc - one value for each restraint(bond distance, angle or dihedral angle sigmas)
virtual void XMLOutput(ostream &os, int indent=0) const
Output to stream in well-formed XML.
vector< MolBond * >::const_iterator FindBond(const MolAtom &, const MolAtom &) const
Searches whether a bond between two atoms already exists.
list< FlipGroup > mvFlipGroup
The list of FlipGroups.
void AddAtom(const REAL x, const REAL y, const REAL z, const ScatteringPower *pPow, const string &name, const bool updateDisplay=true)
Add an atom.
virtual Molecule * CreateCopy() const
vector< MolAtom * >::reverse_iterator FindAtom(const string &name)
Search a MolAtom from its name.
REAL mLogLikelihoodScale
Scale (multiplier) for the log(likelihood)
void BuildStretchModeTwist()
Build the groups of atoms used to twist internally the Molecule, e.g.
virtual const CrystVector_REAL & GetLSQObs(const unsigned int) const
Get the observed values for the LSQ function.
virtual int GetNbComponent() const
Number of components in the scatterer (eg number of point scatterers)
void BuildStretchModeGroups()
Separate StretchMode that break more than their assigned restraint from others.
const RefinableObjClock & GetAtomScattPowClock() const
Get the clock associated to the scattering powers.
const std::vector< MolZAtom > & AsZMatrix(const bool keeporder) const
Molecule as Z-matrix.
list< StretchModeTorsion > mvStretchModeTorsion
List of StretchModeBondLength.
RefinableObjClock mClockRestraint
This clock is the parent of mClockAtomList, mClockBondList, mClockBondAngleList, mClockDihedralAngleL...
void AddRigidGroup(const RigidGroup &, const bool updateDisplay=true)
Add a rigid group of atoms.
RefObjOpt mOptimizeOrientation
Option to optimize the Molecule's orientation.
void OptimizeConformation(const long nbTrial=10000, const REAL stopCost=0.)
Minimize configuration from internal restraints (bond lengths, angles and dihedral angles).
const MolAtom * mpCenterAtom
Atom chosen as center of rotation, if mRotationCenter is set to use an atom rather than the geometric...
virtual void RandomizeConfiguration()
Randomize Configuration (before a global optimization).
void BuildRingList()
Build the list of rings in the molecule.
virtual REAL GetLogLikelihood() const
Get -log(likelihood) of the current configuration for the object.
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.
std::string GetFormula() const
Formula with atoms in alphabetic order.
void RestraintStatus(ostream &os) const
Print the status of all restraints (bond length, angles...)
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.
vector< MolBondAngle * > mvpBondAngle
The list of bond angles.
void RestraintExport(ostream &os) const
Print the restraints (bond length, angles...) as whole labels and number in column text format which ...
void TuneGlobalOptimRotationAmplitude()
Tune the rotation amplitude for free torsions and for the overall Molecule Rotation.
CrystVector_REAL mLSQCalc
Current LSQ Calc - one value for each restraint (bond distance, angle or dihedral angle)
vector< MolBond * >::iterator RemoveBond(const MolBond &, const bool del=true)
Remove a bond.
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.
list< StretchModeBondLength > mvStretchModeBondLength
List of StretchModeBondLength.
virtual string GetComponentName(const int i) const
Name for the i-th component of this scatterer.
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.
void ResetRigidGroupsPar() const
Set the orientation & translation parameters of all rigid groups to 0, after correcting the atomic po...
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.
virtual void Print() const
Print some info about the scatterer (ideally this should be one line...).
std::vector< RigidGroup * >::iterator RemoveRigidGroup(const RigidGroup &group, const bool updateDisplay=true, const bool del=true)
Remove a rigid group of atoms.
const map< MolAtom *, set< MolAtom * > > & GetConnectivityTable()
Get the connectivity table.
void BuildStretchModeBondLength()
Build the groups of atoms moved when stretching a bond length, while respecting the Molecule restrain...
void OptimizeConformationSteepestDescent(const REAL maxStep=0.1, const unsigned nbStep=1)
Optimize the conformation from internal restraints (bond lengths, angles and dihedral angles),...
list< RotorGroup > mvRotorGroupInternal
List of RotorGroups for internal rotations.
void BuildFlipGroup()
Build the groups of atoms that can be flipped.
REAL mLogLikelihood
The current log(likelihood)
virtual const CrystVector_REAL & GetLSQDeriv(const unsigned int n, RefinablePar &par)
Get the first derivative values for the LSQ function, for a given parameter.
Molecule(Crystal &cryst, const string &name="")
Constructor.
vector< MolAtom * > mvpAtom
The list of atoms.
REAL BondAngleRandomChange(const StretchModeBondAngle &mode, const REAL amplitude, const bool respectRestraint=true)
change a bond angle, while respecting the Restraint (if any).
REAL BondLengthRandomChange(const StretchModeBondLength &mode, const REAL amplitude, const bool respectRestraint=true)
Stretch a bond, while respecting the Restraint (if any).
RefinableObjClock & GetBondListClock()
get the clock associated to the list of bonds
virtual void GlobalOptRandomMove(const REAL mutationAmplitude, const RefParType *type)
Make a random move of the current configuration.
RefinableObjClock & GetAtomPositionClock()
Get the clock associated to the atomic positions.
ScatteringComponentList mScattCompList
The list of scattering components.
RefObjOpt mAutoOptimizeConformation
Option to automatically optimize the starting conformation, if the total restraint cost is too high.
std::list< StretchMode * > mvpStretchModeNotFree
Groups of StretchMode breaking restraints (beyond the one they are associated to)
REAL mMDMoveFreq
Frequency of using molecular dynamics move during GlobalOptRandomMove()
virtual const CrystVector_REAL & GetLSQCalc(const unsigned int) const
Get the current calculated value for the LSQ function.
std::set< MolAtom * > mvMDFullAtomGroup
Full list of atoms that can be moved using molecular dynamics This excludes any atom part of a rigid ...
virtual const CrystVector_REAL & GetLSQWeight(const unsigned int) const
Get the weight values for the LSQ function.
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
list< RotorGroup > mvRotorGroupTorsion
List of RotorGroups corresponding to free torsion bonds.
REAL DihedralAngleRandomChange(const StretchModeTorsion &mode, const REAL amplitude, const bool respectRestraint=true)
Change a dihedral angle, while respecting the Restraint (if any).
Defines a group of atoms which can be rotated around an axis defined by two other atoms.
RotorGroup(const MolAtom &at1, const MolAtom &at2)
Constructor, with the two atoms around which the rotation shall be made.
When 3(A1..1n) or more atoms are connected to a same atom A, it defines a 'flip' group,...
list< pair< const MolAtom *, set< MolAtom * > > > mvRotatedChainList
The set of atoms that are to be rotated during the flip.
const MolAtom * mpAtom2
The second atom defining the rotation axis.
const MolAtom * mpAtom0
The atom which is an asymmetric center.
FlipGroup(const MolAtom &at0, const MolAtom &at1, const MolAtom &at2)
Constructor, with the central atom.
const MolAtom * mpAtom1
The first atom defining the rotation axis.
Crystal * mpCryst
The crystal in which the Scatterer is This is needed so that we can know which scattering powers are ...
virtual void SetZ(const REAL z)
Z coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
virtual void SetY(const REAL y)
Y coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
REAL GetOccupancy() const
Get the occupancy of the scatterer (0.
const Crystal & GetCrystal() const
In which crystal is this Scatterer included ?
RefinableObjClock mClockScattCompList
CrystVector_REAL mXYZ
coordinates of the scatterer (or of its center..)
REAL GetX() const
X coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
REAL GetZ() const
Z coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
RefinableObjClock mClockScatterer
Last time anything (number of atoms, positions, scattering power) was changed.
REAL mOccupancy
Occupancy : 0 <= occ <= 1 For a multi-atom scatterer (polyhedron,..), this is the overall occupancy o...
REAL GetY() const
Y coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
virtual void SetX(const REAL x)
X coordinate (fractionnal) of the scatterer (for complex scatterers, this corresponds to the position...
Abstract Base Class to describe the scattering power of any Scatterer component in a crystal.
virtual REAL GetRadius() const =0
Return the physical radius of this type of scatterer (for 3D display purposes).
The Scattering Power for an Atom.
int GetAtomicNumber() const
Atomic number for this atom.
list of scattering positions in a crystal, associated with the corresponding occupancy and a pointer ...
long GetNbComponent() const
Number of components.
CrystVector_REAL GetLatticePar() const
Lattice parameters (a,b,c,alpha,beta,gamma) as a 6-element vector in Angstroems and radians.
void OrthonormalToFractionalCoords(REAL &x, REAL &y, REAL &z) const
Get fractional cartesian coordinates for a set of (x,y,z) orthonormal coordinates.
void FractionalToOrthonormalCoords(REAL &x, REAL &y, REAL &z) const
Get orthonormal cartesian coordinates for a set of (x,y,z) fractional coordinates.
ZScatterer: the basic type of complex scatterers, where atom positions are defined using a standard "...
REAL GetZAtomX(const int i) const
Get the X fractionnal coordinate of atom i.
const ObjRegistry< ZAtom > & GetZAtomRegistry() const
Access to the registry of ZAtoms.
long GetZBondAtom(const int i) const
Index of the 1st atom used to define the i-th atom in the Z-Matrix (the one from which the bondlength...
REAL GetZAtomZ(const int i) const
Get the Z fractionnal coordinate of atom i.
long GetZDihedralAngleAtom(const int i) const
Index of the 3rd atom used to define the i-th atom in the Z-Matrix (the one from which the dihedral a...
REAL GetZAtomY(const int i) const
Get the Y fractionnal coordinate of atom i.
long GetZAngleAtom(const int i) const
Index of the 2nd atom used to define the i-th atom in the Z-Matrix (the one from which the angle is c...
Vector library (Blitz++ mimic) for ObjCryst++.
Cubic spline interpolation.
void AddRefinableObj(RefinableObj &)
Add a refined object. All sub-objects are also added.
Base object for Monte-Carlo Global Optimization methods.
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.
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.
class to input or output a well-formatted xml beginning or ending tag.
class of refinable parameter types.
bool IsDescendantFromOrSameAs(const RefParType *type) const
Returns true if the parameter is a descendant of 'type'.
We need to record exactly when refinable objects have been modified for the last time (to avoid re-co...
void Reset()
Reset a Clock to 0, to force an update.
void AddChild(const RefinableObjClock &)
Add a 'child' clock.
void Click()
Record an event for this clock (generally, the 'time' an object has been modified,...
Generic class for parameters of refinable objects.
void XMLOutput(ostream &os, const string &name, int indent=0) const
XMLOutput to stream in well-formed XML.
void SetDerivStep(const REAL)
Fixed step to use to compute numerical derivative.
REAL GetMin() const
Minimum value allowed (if limited or periodic)
REAL GetValue() const
of the parameter.
void SetGlobalOptimStep(const REAL)
Maximum step to use during Global Optimization algorithms.
void AssignClock(RefinableObjClock &clock)
REAL GetMax() const
Get the maximum value allowed (if limited)
void SetName(const string &)
Set the name of the parameter. It should be unique in the RefinableObj.
void XMLInput(istream &is, const XMLCrystTag &tag)
XMLInput From stream.
void XMLOutput(ostream &os, int indent=0) const
XMLOutput to stream in well-formed XML.
virtual void EndOptimization()
This should be called by any optimization class at the end of an optimization.
void RestoreParamSet(const unsigned long id)
Restore a saved set of values.
void AddPar(const RefinablePar &newRefPar)
Add a refinable parameter.
void AddRestraint(Restraint *pNewRestraint)
Add a new restraint.
bool IsBeingRefined() const
Is the object being refined ? (Can be refined by one algorithm at a time only.)
virtual void SetName(const string &name)
Name of the object.
virtual const CrystVector_REAL & GetLSQDeriv(const unsigned int, RefinablePar &)
Get the first derivative values for the LSQ function, for a given parameter.
RefinablePar & GetPar(const long i)
Access all parameters in the order they were inputted.
vector< Restraint * >::iterator RemoveRestraint(Restraint *pRestraint)
Remove a restraint from the list of known restraints.
vector< RefinablePar * >::iterator RemovePar(RefinablePar *refPar)
Remove a refinable parameter.
unsigned int GetNbOption() const
Number of Options for this object.
virtual const string & GetName() const
Name of the object.
void SaveParamSet(const unsigned long id) const
Save the current set of refined values over a previously-created set of saved values.
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.
RefObjOpt & GetOption(const unsigned int i)
Access to the options.
ObjRegistry< RefObjOpt > mOptionRegistry
List of options for this object.
void AddOption(RefObjOpt *opt)
void ClearParamSet(const unsigned long id) const
Erase the param set with the given id, releasing memory.
unsigned long CreateParamSet(const string name="") const
Save the current set of refined values in a new set.
virtual void UpdateDisplay() const
If there is an interface, this should be automatically be called each time there is a 'new,...
virtual REAL GetLogLikelihood() const
Get -log(likelihood) of the current configuration for the object.
vector< Restraint * > mvpRestraint
Vector of pointers to the restraints for this object.
string mName
Name for this RefinableObject. Should be unique, at least in the same scope.+.
virtual void RandomizeConfiguration()
Randomize Configuration (before a global optimization).
virtual void GlobalOptRandomMove(const REAL mutationAmplitude, const RefParType *type=gpRefParTypeObjCryst)
Make a random move of the current configuration.
Abstract base class for all objects in wxCryst.
wx class for MolAtom objects
wx class for MolBond objects
wx class for MolBondAngle objects
wx class for MolDihedralAngle objects
wxCryst class for Molecule objects