30 void DRMSD::setFromPDB(
const PDB&pdb,
double lbound,
double ubound){
38 void DRMSD::setReference(
const vector<Vector> & reference,
double lbound,
double ubound)
40 natoms = reference.size();
41 for(
unsigned i=0;i<natoms-1;++i){
42 for(
unsigned j=i+1;j<natoms;++j){
43 double distance =
delta(reference[i],reference[j]).modulo();
44 if(distance > lbound && distance < ubound){
45 targets[make_pair(i,j)] = distance;
51 double DRMSD::calculate(
const std::vector<Vector> & positions,
52 std::vector<Vector> &derivatives,
Tensor& virial)
const {
55 return DRMSD::calculate(positions,fake_pbc,derivatives,virial,
false);
58 double DRMSD::calculate(
const std::vector<Vector> & positions,
const Pbc& pbc,
59 std::vector<Vector> &derivatives,
Tensor& virial,
bool do_pbc)
const {
61 plumed_assert(positions.size()==natoms && derivatives.size()==natoms );
64 double drmsd=0.; virial.
zero();
65 for(
unsigned i=0;i<derivatives.size();++i) derivatives[i].zero();
66 for(map< pair <unsigned,unsigned> ,
double>::const_iterator it=targets.begin();it!=targets.end();++it){
67 unsigned i=it->first.first;
68 unsigned j=it->first.second;
69 if(do_pbc){distance=pbc.
distance( positions[i] , positions[j] );}
70 else{distance=
delta( positions[i] , positions[j] );}
71 double len = distance.modulo();
72 double diff = len - it->second;
74 derivatives[i] += -( diff / len ) * distance;
75 derivatives[j] += ( diff / len ) * distance;
76 virial -= ( diff / len ) *
Tensor(distance,distance);
79 double npairs =
static_cast<double>(targets.size());
81 drmsd=sqrt( drmsd / npairs );
83 double idrmsd=1.0/( drmsd * npairs );
87 for(
unsigned i=0;i<natoms;++i){derivatives[i] *= idrmsd;}
void zero()
set it to zero
Class implementing fixed size matrices of doubles.
Class implementing fixed size vectors of doubles.
Vector distance(const Vector &, const Vector &, int *nshifts) const
internal version of distance, also returns the number of attempted shifts (used in Pbc::test())...
const std::vector< Vector > & getPositions() const
Access to the position array.
VectorGeneric< n > delta(const VectorGeneric< n > &v1, const VectorGeneric< n > &v2)