Skip to content

Commit

Permalink
New interface for calculateCV to make it easier to change the input a…
Browse files Browse the repository at this point in the history 
…rguments in future

We now pass an object that contains all the input arguments to calculateCV.  This means
that if we need to pass new things to this command in future you wont need to modify as much code
  • Loading branch information
Gareth Aneurin Tribello authored and Gareth Aneurin Tribello committed Feb 5, 2025
1 parent 6fa99dc commit 2986098
Show file tree
Hide file tree
Showing 11 changed files with 123 additions and 95 deletions.
16 changes: 7 additions & 9 deletions src/colvar/Angle.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ Calculate multiple angles.

class Angle : public Colvar {
bool pbc;
std::vector<double> value, masses, charges;
std::vector<double> value;
std::vector<std::vector<Vector> > derivs;
std::vector<Tensor> virial;
public:
Expand All @@ -110,8 +110,7 @@ class Angle : public Colvar {
static void registerKeywords( Keywords& keys );
static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
static unsigned getModeAndSetupValues( ActionWithValue* av );
static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa );
};

Expand Down Expand Up @@ -169,23 +168,22 @@ Angle::Angle(const ActionOptions&ao):
void Angle::calculate() {

if(pbc) makeWhole();
calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
setValue( value[0] );
for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
setBoxDerivatives( virial[0] );
}

void Angle::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
void Angle::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
Vector dij,dik;
dij=delta(pos[2],pos[3]);
dik=delta(pos[1],pos[0]);
dij=delta(cvin.pos[2],cvin.pos[3]);
dik=delta(cvin.pos[1],cvin.pos[0]);
Vector ddij,ddik; PLMD::Angle a;
vals[0]=a.compute(dij,dik,ddij,ddik);
derivs[0][0]=ddik; derivs[0][1]=-ddik;
derivs[0][2]=-ddij; derivs[0][3]=ddij;
setBoxDerivativesNoPbc( pos, derivs, virial );
setBoxDerivativesNoPbc( cvin.pos, derivs, virial );
}

}
Expand Down
22 changes: 10 additions & 12 deletions src/colvar/DihedralCorrelation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ Measure the correlation between a multiple pairs of dihedral angles
class DihedralCorrelation : public Colvar {
private:
bool pbc;
std::vector<double> value, masses, charges;
std::vector<double> value;
std::vector<std::vector<Vector> > derivs;
std::vector<Tensor> virial;
public:
Expand All @@ -73,8 +73,7 @@ class DihedralCorrelation : public Colvar {
static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
static unsigned getModeAndSetupValues( ActionWithValue* av );
void calculate() override;
static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa );
};

Expand Down Expand Up @@ -126,26 +125,25 @@ unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
void DihedralCorrelation::calculate() {

if(pbc) makeWhole();
calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
setValue( value[0] );
for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
setBoxDerivatives( virial[0] );
}

void DihedralCorrelation::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
void DihedralCorrelation::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
const Vector d10=delta(pos[1],pos[0]);
const Vector d11=delta(pos[2],pos[1]);
const Vector d12=delta(pos[3],pos[2]);
const Vector d10=delta(cvin.pos[1],cvin.pos[0]);
const Vector d11=delta(cvin.pos[2],cvin.pos[1]);
const Vector d12=delta(cvin.pos[3],cvin.pos[2]);

Vector dd10,dd11,dd12;
PLMD::Torsion t1;
const double phi1 = t1.compute(d10,d11,d12,dd10,dd11,dd12);

const Vector d20=delta(pos[5],pos[4]);
const Vector d21=delta(pos[6],pos[5]);
const Vector d22=delta(pos[7],pos[6]);
const Vector d20=delta(cvin.pos[5],cvin.pos[4]);
const Vector d21=delta(cvin.pos[6],cvin.pos[5]);
const Vector d22=delta(cvin.pos[7],cvin.pos[6]);

Vector dd20,dd21,dd22;
PLMD::Torsion t2;
Expand Down
37 changes: 17 additions & 20 deletions src/colvar/Dipole.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ class Dipole : public Colvar {
std::vector<AtomNumber> ga_lista;
bool components;
bool nopbc;
std::vector<double> value, masses, charges;
std::vector<double> value;
std::vector<std::vector<Vector> > derivs;
std::vector<Tensor> virial;
Value* valuex=nullptr;
Expand All @@ -94,8 +94,7 @@ class Dipole : public Colvar {
static unsigned getModeAndSetupValues( ActionWithValue* av );
void calculate() override;
static void registerKeywords(Keywords& keys);
static void calculateCV( const unsigned& mode, const std::vector<double>& masses, std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa );
};

Expand Down Expand Up @@ -123,7 +122,7 @@ Dipole::Dipole(const ActionOptions&ao):
derivs(1),
virial(1)
{
parseAtomList(-1,ga_lista,this); charges.resize(ga_lista.size());
parseAtomList(-1,ga_lista,this);
components=(getModeAndSetupValues(this)==1);
if( components ) {
value.resize(3); derivs.resize(3); virial.resize(3);
Expand Down Expand Up @@ -166,17 +165,18 @@ unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) {
// calculator
void Dipole::calculate()
{
if( !chargesWereSet ) error("charges were not set by MD code");

if(!nopbc) makeWhole();
unsigned N=getNumberOfAtoms();
for(unsigned i=0; i<N; ++i) charges[i]=getCharge(i);

if(!components) {
calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
for(unsigned i=0; i<N; i++) setAtomsDerivatives(i,derivs[0][i]);
setBoxDerivatives(virial[0]);
setValue(value[0]);
} else {
calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), value, derivs, virial, this );
for(unsigned i=0; i<N; i++) {
setAtomsDerivatives(valuex,i,derivs[0][i]);
setAtomsDerivatives(valuey,i,derivs[1][i]);
Expand All @@ -191,34 +191,31 @@ void Dipole::calculate()
}
}

void Dipole::calculateCV( const unsigned& mode, const std::vector<double>& masses, std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
void Dipole::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
unsigned N=pos.size(); double ctot=0.;
for(unsigned i=0; i<N; ++i) ctot += charges[i];
unsigned N=cvin.pos.size(); double ctot=0.;
for(unsigned i=0; i<N; ++i) ctot += cvin.charges[i];
ctot/=(double)N;

Vector dipje;
for(unsigned i=0; i<N; ++i) {
charges[i]-=ctot; dipje += charges[i]*pos[i];
}
for(unsigned i=0; i<N; ++i) dipje += (cvin.charges[i]-ctot)*cvin.pos[i];

if( mode==1 ) {
if( cvin.mode==1 ) {
for(unsigned i=0; i<N; i++) {
derivs[0][i]=charges[i]*Vector(1.0,0.0,0.0);
derivs[1][i]=charges[i]*Vector(0.0,1.0,0.0);
derivs[2][i]=charges[i]*Vector(0.0,0.0,1.0);
derivs[0][i]=(cvin.charges[i]-ctot)*Vector(1.0,0.0,0.0);
derivs[1][i]=(cvin.charges[i]-ctot)*Vector(0.0,1.0,0.0);
derivs[2][i]=(cvin.charges[i]-ctot)*Vector(0.0,0.0,1.0);
}
for(unsigned i=0; i<3; ++i ) vals[i] = dipje[i];
} else {
vals[0] = dipje.modulo();
double idip = 1./vals[0];
for(unsigned i=0; i<N; i++) {
double dfunc=charges[i]*idip;
double dfunc=(cvin.charges[i]-ctot)*idip;
derivs[0][i] = dfunc*dipje;
}
}
setBoxDerivativesNoPbc( pos, derivs, virial );
setBoxDerivativesNoPbc( cvin.pos, derivs, virial );
}

}
Expand Down
24 changes: 11 additions & 13 deletions src/colvar/Distance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,7 @@ class Distance : public Colvar {
bool scaled_components;
bool pbc;

std::vector<double> value, masses, charges;
std::vector<double> value;
std::vector<std::vector<Vector> > derivs;
std::vector<Tensor> virial;
public:
Expand All @@ -138,8 +138,7 @@ class Distance : public Colvar {
static unsigned getModeAndSetupValues( ActionWithValue* av );
// active methods:
void calculate() override;
static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa );
};

Expand Down Expand Up @@ -227,7 +226,7 @@ void Distance::calculate() {
if(pbc) makeWhole();

if( components ) {
calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), value, derivs, virial, this );
Value* valuex=getPntrToComponent("x");
Value* valuey=getPntrToComponent("y");
Value* valuez=getPntrToComponent("z");
Expand All @@ -244,7 +243,7 @@ void Distance::calculate() {
setBoxDerivatives(valuez,virial[2]);
valuez->set(value[2]);
} else if( scaled_components ) {
calculateCV( 2, masses, charges, getPositions(), value, derivs, virial, this );
calculateCV( ColvarInput::createColvarInput( 2, getPositions(), this ), value, derivs, virial, this );

Value* valuea=getPntrToComponent("a");
Value* valueb=getPntrToComponent("b");
Expand All @@ -256,21 +255,20 @@ void Distance::calculate() {
for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuec,i,derivs[2][i] );
valuec->set(value[2]);
} else {
calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
for(unsigned i=0; i<2; ++i) setAtomsDerivatives(i,derivs[0][i] );
setBoxDerivatives(virial[0]);
setValue (value[0]);
}
}

void Distance::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
void Distance::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
Vector distance=delta(pos[0],pos[1]);
Vector distance=delta(cvin.pos[0],cvin.pos[1]);
const double value=distance.modulo();
const double invvalue=1.0/value;

if(mode==1) {
if(cvin.mode==1) {
derivs[0][0] = Vector(-1,0,0);
derivs[0][1] = Vector(+1,0,0);
vals[0] = distance[0];
Expand All @@ -282,8 +280,8 @@ void Distance::calculateCV( const unsigned& mode, const std::vector<double>& mas
derivs[2][0] = Vector(0,0,-1);
derivs[2][1] = Vector(0,0,+1);
vals[2] = distance[2];
setBoxDerivativesNoPbc( pos, derivs, virial );
} else if(mode==2) {
setBoxDerivativesNoPbc( cvin.pos, derivs, virial );
} else if(cvin.mode==2) {
Vector d=aa->getPbc().realToScaled(distance);
derivs[0][0] = matmul(aa->getPbc().getInvBox(),Vector(-1,0,0));
derivs[0][1] = matmul(aa->getPbc().getInvBox(),Vector(+1,0,0));
Expand All @@ -297,7 +295,7 @@ void Distance::calculateCV( const unsigned& mode, const std::vector<double>& mas
} else {
derivs[0][0] = -invvalue*distance;
derivs[0][1] = invvalue*distance;
setBoxDerivativesNoPbc( pos, derivs, virial );
setBoxDerivativesNoPbc( cvin.pos, derivs, virial );
vals[0] = value;
}
}
Expand Down
17 changes: 17 additions & 0 deletions src/colvar/MultiColvarTemplate.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,5 +20,22 @@
along with plumed. If not, see <http://www.gnu.org/licenses/>.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "MultiColvarTemplate.h"
#include "core/Colvar.h"

namespace PLMD {
namespace colvar {

ColvarInput::ColvarInput( const unsigned& m, const std::vector<Vector>& p, const std::vector<double>& w, const std::vector<double>& q ) :
mode(m),
pos(p),
mass(w),
charges(q)
{
}

ColvarInput ColvarInput::createColvarInput( const unsigned& m, const std::vector<Vector>& p, const Colvar* colv ) {
return ColvarInput( m, p, colv->getMasses(), colv->getCharges(true) );
}

}
}
15 changes: 14 additions & 1 deletion src/colvar/MultiColvarTemplate.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,21 @@
#include "core/ParallelTaskManager.h"

namespace PLMD {

class Colvar;

namespace colvar {

class ColvarInput {
public:
unsigned mode;
const std::vector<Vector>& pos;
const std::vector<double>& mass;
const std::vector<double>& charges;
ColvarInput( const unsigned& m, const std::vector<Vector>& p, const std::vector<double>& w, const std::vector<double>& q );
static ColvarInput createColvarInput( const unsigned& m, const std::vector<Vector>& p, const Colvar* colv );
};

template <class T>
class MultiColvarTemplate : public ActionWithVector {
private:
Expand Down Expand Up @@ -164,7 +177,7 @@ void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue
std::vector<Tensor> & virial( myvals.getFirstAtomVirialVector() );
std::vector<std::vector<Vector> > & derivs( myvals.getFirstAtomDerivativeVector() );
// Calculate the CVs using the method in the Colvar
T::calculateCV( mode, mass, charge, fpositions, values, derivs, virial, this );
T::calculateCV( ColvarInput(mode, fpositions, mass, charge), values, derivs, virial, this );
for(unsigned i=0; i<values.size(); ++i) myvals.setValue( i, values[i] );
// Finish if there are no derivatives
if( doNotCalculateDerivatives() ) return;
Expand Down
14 changes: 6 additions & 8 deletions src/colvar/Plane.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ namespace colvar {
class Plane : public Colvar {
private:
bool pbc;
std::vector<double> value, masses, charges;
std::vector<double> value;
std::vector<std::vector<Vector> > derivs;
std::vector<Tensor> virial;
public:
Expand All @@ -72,8 +72,7 @@ class Plane : public Colvar {
static unsigned getModeAndSetupValues( ActionWithValue* av );
// active methods:
void calculate() override;
static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa );
};

Expand Down Expand Up @@ -135,17 +134,16 @@ Plane::Plane(const ActionOptions&ao):
void Plane::calculate() {

if(pbc) makeWhole();
calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
setValue( value[0] );
for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
setBoxDerivatives( virial[0] );
}

void Plane::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
void Plane::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
Vector d1=delta( pos[1], pos[0] );
Vector d2=delta( pos[2], pos[3] );
Vector d1=delta( cvin.pos[1], cvin.pos[0] );
Vector d2=delta( cvin.pos[2], cvin.pos[3] );
Vector cp = crossProduct( d1, d2 );

derivs[0][0] = crossProduct( Vector(-1.0,0,0), d2 );
Expand Down
Loading

0 comments on commit 2986098

Please sign in to comment.