From 8f1accede393dd17caf83c0e535218554db71524 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Tue, 19 Nov 2024 11:28:30 +0100 Subject: [PATCH 01/10] Macro for fasterMulticolvar+named modes --- src/colvar/Angle.cpp | 24 ++++--- src/colvar/DihedralCorrelation.cpp | 19 +++-- src/colvar/Dipole.cpp | 41 +++++++---- src/colvar/Distance.cpp | 51 ++++++++------ src/colvar/MultiColvarTemplate.h | 109 +++++++++++++++++++++-------- src/colvar/Plane.cpp | 16 ++--- src/colvar/Position.cpp | 28 ++++---- src/colvar/SelectMassCharge.cpp | 22 +++--- src/colvar/Torsion.cpp | 52 +++++++++----- src/core/ActionWithVector.cpp | 10 ++- src/core/Value.cpp | 4 +- src/crystdistrib/Quaternion.cpp | 16 ++--- 12 files changed, 239 insertions(+), 153 deletions(-) diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp index 846e0432f7..0a6579cd00 100644 --- a/src/colvar/Angle.cpp +++ b/src/colvar/Angle.cpp @@ -104,15 +104,13 @@ class Angle : public Colvar { std::vector > derivs; std::vector virial; public: + MULTICOLVAR_SETTINGS(multiColvars::emptyMode); explicit Angle(const ActionOptions&); // active methods: void calculate() override; static void registerKeywords( Keywords& keys ); static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); - static void calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); + }; typedef ColvarShortcut AngleShortcut; @@ -140,8 +138,10 @@ void Angle::parseAtomList( const int& num, std::vector& atoms, Actio } else if( num<0 || atoms.size()>0 ) aa->error("Number of specified atoms should be either 3 or 4"); } -unsigned Angle::getModeAndSetupValues( ActionWithValue* av ) { - av->addValueWithDerivatives(); av->setNotPeriodic(); return 0; +Angle::Modetype Angle::getModeAndSetupValues( ActionWithValue* av ) { + av->addValueWithDerivatives(); + av->setNotPeriodic(); + return {}; } Angle::Angle(const ActionOptions&ao): @@ -169,19 +169,21 @@ Angle::Angle(const ActionOptions&ao): void Angle::calculate() { if(pbc) makeWhole(); - calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this ); + calculateCV( {}, masses, charges, getPositions(), value, derivs, virial, this ); setValue( value[0] ); - for(unsigned i=0; i& masses, const std::vector& charges, +void Angle::calculateCV( Modetype /*mode*/, const std::vector& /*masses*/, const std::vector& /*charges*/, const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { + std::vector& virial, const ActionAtomistic* /*aa*/ ) { Vector dij,dik; dij=delta(pos[2],pos[3]); dik=delta(pos[1],pos[0]); - Vector ddij,ddik; PLMD::Angle a; + 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; diff --git a/src/colvar/DihedralCorrelation.cpp b/src/colvar/DihedralCorrelation.cpp index 19e08feb1b..28ff787af9 100644 --- a/src/colvar/DihedralCorrelation.cpp +++ b/src/colvar/DihedralCorrelation.cpp @@ -68,14 +68,11 @@ class DihedralCorrelation : public Colvar { std::vector > derivs; std::vector virial; public: + MULTICOLVAR_SETTINGS(multiColvars::emptyMode); static void registerKeywords( Keywords& keys ); explicit DihedralCorrelation(const ActionOptions&); static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); void calculate() override; - static void calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); }; typedef ColvarShortcut DihedralCorrelationShortcut; @@ -119,22 +116,24 @@ void DihedralCorrelation::parseAtomList( const int& num, std::vector } } -unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) { - av->addValueWithDerivatives(); av->setNotPeriodic(); return 0; +DihedralCorrelation::Modetype DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) { + av->addValueWithDerivatives(); + av->setNotPeriodic(); + return {}; } void DihedralCorrelation::calculate() { if(pbc) makeWhole(); - calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this ); + calculateCV( {}, masses, charges, getPositions(), value, derivs, virial, this ); setValue( value[0] ); for(unsigned i=0; i& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { +void DihedralCorrelation::calculateCV(Modetype /*mode*/, const std::vector& masses, const std::vector& charges, + const std::vector& pos, std::vector& vals, std::vector >& derivs, + std::vector& 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]); diff --git a/src/colvar/Dipole.cpp b/src/colvar/Dipole.cpp index e0a9212455..10873ca79c 100644 --- a/src/colvar/Dipole.cpp +++ b/src/colvar/Dipole.cpp @@ -91,12 +91,19 @@ class Dipole : public Colvar { public: explicit Dipole(const ActionOptions&); static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); void calculate() override; static void registerKeywords(Keywords& keys); - static void calculateCV( const unsigned& mode, const std::vector& masses, std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); + MULTICOLVAR_SETTINGS_MODE(multiColvars::components); + MULTICOLVAR_SETTINGS_SETUPF(); + //Declaring this playinly becasue we are using modifying the charges + static void calculateCV( Modetype mode, + const std::vector& masses, + std::vector& charges, + const std::vector& pos, + std::vector& vals, + std::vector >& derivs, + std::vector& virial, + const ActionAtomistic* aa ); }; typedef ColvarShortcut DipoleShortcut; @@ -123,8 +130,9 @@ Dipole::Dipole(const ActionOptions&ao): derivs(1), virial(1) { - parseAtomList(-1,ga_lista,this); charges.resize(ga_lista.size()); - components=(getModeAndSetupValues(this)==1); + parseAtomList(-1,ga_lista,this); + charges.resize(ga_lista.size()); + components=(getModeAndSetupValues(this)==Modetype::withCompontents); if( components ) { value.resize(3); derivs.resize(3); virial.resize(3); valuex=getPntrToComponent("x"); @@ -152,15 +160,17 @@ void Dipole::parseAtomList( const int& num, std::vector& t, ActionAt } } -unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) { - bool c; av->parseFlag("COMPONENTS",c); +Dipole::Modetype Dipole::getModeAndSetupValues( ActionWithValue* av ) { + bool c; + av->parseFlag("COMPONENTS",c); if( c ) { av->addComponentWithDerivatives("x"); av->componentIsNotPeriodic("x"); av->addComponentWithDerivatives("y"); av->componentIsNotPeriodic("y"); av->addComponentWithDerivatives("z"); av->componentIsNotPeriodic("z"); - return 1; + return Modetype::withCompontents; } - av->addValueWithDerivatives(); av->setNotPeriodic(); return 0; + av->addValueWithDerivatives(); av->setNotPeriodic(); + return Modetype::noComponents; } // calculator @@ -171,12 +181,12 @@ void Dipole::calculate() for(unsigned i=0; i& masses, std::vector& charges, +void Dipole::calculateCV( Modetype mode, const std::vector& masses, std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) { unsigned N=pos.size(); double ctot=0.; @@ -200,10 +210,11 @@ void Dipole::calculateCV( const unsigned& mode, const std::vector& masse Vector dipje; for(unsigned i=0; i& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); // active methods: void calculate() override; - static void calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); + MULTICOLVAR_SETTINGS(multiColvars::scaledComponents); }; typedef ColvarShortcut DistanceShortcut; @@ -186,11 +183,18 @@ Distance::Distance(const ActionOptions&ao): if(pbc) log.printf(" using periodic boundary conditions\n"); else log.printf(" without periodic boundary conditions\n"); - unsigned mode = getModeAndSetupValues( this ); - if(mode==1) components=true; else if(mode==2) scaled_components=true; + Modetype mode = getModeAndSetupValues( this ); + if(mode==Modetype::withCompontents) { + components=true; + } else if(mode==Modetype::scaledComponents) { + scaled_components=true; + } if( components || scaled_components ) { - value.resize(3); derivs.resize(3); virial.resize(3); - for(unsigned i=0; i<3; ++i) derivs[i].resize(2); + value.resize(3); + derivs.resize(3); + virial.resize(3); + for(unsigned i=0; i<3; ++i) + derivs[i].resize(2); } requestAtoms(atoms); } @@ -200,25 +204,28 @@ void Distance::parseAtomList( const int& num, std::vector& t, Action if( t.size()==2 ) aa->log.printf(" between atoms %d %d\n",t[0].serial(),t[1].serial()); } -unsigned Distance::getModeAndSetupValues( ActionWithValue* av ) { - bool c; av->parseFlag("COMPONENTS",c); - bool sc; av->parseFlag("SCALED_COMPONENTS",sc); - if( c && sc ) av->error("COMPONENTS and SCALED_COMPONENTS are not compatible"); - +Distance::Modetype Distance::getModeAndSetupValues( ActionWithValue* av ) { + bool c; + av->parseFlag("COMPONENTS",c); + bool sc; + av->parseFlag("SCALED_COMPONENTS",sc); + if( c && sc ) { + av->error("COMPONENTS and SCALED_COMPONENTS are not compatible"); + } if(c) { av->addComponentWithDerivatives("x"); av->componentIsNotPeriodic("x"); av->addComponentWithDerivatives("y"); av->componentIsNotPeriodic("y"); av->addComponentWithDerivatives("z"); av->componentIsNotPeriodic("z"); av->log<<" WARNING: components will not have the proper periodicity - see manual\n"; - return 1; + return Modetype::withCompontents; } else if(sc) { av->addComponentWithDerivatives("a"); av->componentIsPeriodic("a","-0.5","+0.5"); av->addComponentWithDerivatives("b"); av->componentIsPeriodic("b","-0.5","+0.5"); av->addComponentWithDerivatives("c"); av->componentIsPeriodic("c","-0.5","+0.5"); - return 2; + return Modetype::scaledComponents; } av->addValueWithDerivatives(); av->setNotPeriodic(); - return 0; + return Modetype::noComponents; } // calculator @@ -227,7 +234,7 @@ void Distance::calculate() { if(pbc) makeWhole(); if( components ) { - calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this ); + calculateCV( Modetype::withCompontents, masses, charges, getPositions(), value, derivs, virial, this ); Value* valuex=getPntrToComponent("x"); Value* valuey=getPntrToComponent("y"); Value* valuez=getPntrToComponent("z"); @@ -244,7 +251,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( Modetype::scaledComponents, masses, charges, getPositions(), value, derivs, virial, this ); Value* valuea=getPntrToComponent("a"); Value* valueb=getPntrToComponent("b"); @@ -256,21 +263,21 @@ 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( Modetype::noComponents, masses, charges, getPositions(), 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& masses, const std::vector& charges, +void Distance::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) { Vector distance=delta(pos[0],pos[1]); const double value=distance.modulo(); const double invvalue=1.0/value; - if(mode==1) { + if(mode==Modetype::withCompontents) { derivs[0][0] = Vector(-1,0,0); derivs[0][1] = Vector(+1,0,0); vals[0] = distance[0]; @@ -283,7 +290,7 @@ void Distance::calculateCV( const unsigned& mode, const std::vector& mas derivs[2][1] = Vector(0,0,+1); vals[2] = distance[2]; setBoxDerivativesNoPbc( pos, derivs, virial ); - } else if(mode==2) { + } else if(mode==Modetype::scaledComponents) { 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)); diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h index bd11210de5..578aedcbc7 100644 --- a/src/colvar/MultiColvarTemplate.h +++ b/src/colvar/MultiColvarTemplate.h @@ -27,11 +27,47 @@ namespace PLMD { namespace colvar { -template +namespace multiColvars { +struct emptyMode {}; +enum class components { + withCompontents, + noComponents +}; +enum class plainOrScaled { + scaled, + plain +}; +enum class scaledComponents { + withCompontents, + scaledComponents, + noComponents +}; +} +#define MULTICOLVAR_SETTINGS_MODE(type) using Modetype=type +#define MULTICOLVAR_SETTINGS_SETUPF() static Modetype getModeAndSetupValues( ActionWithValue* av ) +#define MULTICOLVAR_SETTINGS_CALCULATE_CONST() static void calculateCV( Modetype mode, \ + const std::vector& masses, \ + const std::vector& charges, \ + const std::vector& pos, \ + std::vector& vals, \ + std::vector >& derivs,\ + std::vector& virial, \ + const ActionAtomistic* aa ) + +#define MULTICOLVAR_SETTINGS(type) MULTICOLVAR_SETTINGS_MODE(type); \ + MULTICOLVAR_SETTINGS_SETUPF(); \ + MULTICOLVAR_SETTINGS_CALCULATE_CONST() + + + +//^no ';' here so that the macro will "consume" it and not generate a double semicolon error + +template class MultiColvarTemplate : public ActionWithVector { private: + using Modetype =typename Colvar::Modetype; /// An index that decides what we are calculating - unsigned mode; + Modetype mode; /// Are we using pbc to calculate the CVs bool usepbc; /// Do we reassemble the molecule @@ -48,9 +84,9 @@ class MultiColvarTemplate : public ActionWithVector { void calculate() override; }; -template -void MultiColvarTemplate::registerKeywords(Keywords& keys ) { - T::registerKeywords( keys ); +template +void MultiColvarTemplate::registerKeywords(Keywords& keys ) { + Colvar::registerKeywords( keys ); keys.add("optional","MASK","the label for a sparse matrix that should be used to determine which elements of the matrix should be computed"); unsigned nkeys = keys.size(); for(unsigned i=0; i::registerKeywords(Keywords& keys ) { if( keys.outputComponentExists(".#!value") ) keys.setValueDescription("the " + keys.getDisplayName() + " for each set of specified atoms"); } -template -MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): +template +MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): Action(ao), ActionWithVector(ao), - mode(0), + mode(), usepbc(true), wholemolecules(false) { @@ -75,7 +111,7 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): } else { std::vector t; for(int i=1;; ++i ) { - T::parseAtomList( i, t, this ); + Colvar::parseAtomList( i, t, this ); if( t.empty() ) break; if( i==1 ) { ablocks.resize(t.size()); } @@ -103,36 +139,40 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): else log.printf(" without periodic boundary conditions\n"); // Setup the values - mode = T::getModeAndSetupValues( this ); + mode = Colvar::getModeAndSetupValues( this ); } -template -unsigned MultiColvarTemplate::getNumberOfDerivatives() { +template +unsigned MultiColvarTemplate::getNumberOfDerivatives() { return 3*getNumberOfAtoms()+9; } -template -void MultiColvarTemplate::calculate() { +template +void MultiColvarTemplate::calculate() { if( wholemolecules ) makeWhole(); runAllTasks(); } -template -void MultiColvarTemplate::addValueWithDerivatives( const std::vector& shape ) { +template +void MultiColvarTemplate::addValueWithDerivatives( const std::vector& shape ) { std::vector s(1); s[0]=ablocks[0].size(); addValue( s ); } -template -void MultiColvarTemplate::addComponentWithDerivatives( const std::string& name, const std::vector& shape ) { +template +void MultiColvarTemplate::addComponentWithDerivatives( const std::string& name, const std::vector& shape ) { std::vector s(1); s[0]=ablocks[0].size(); addComponent( name, s ); } -template -void MultiColvarTemplate::performTask( const unsigned& task_index, MultiValue& myvals ) const { +template +void MultiColvarTemplate::performTask( const unsigned& task_index, MultiValue& myvals ) const { // Retrieve the positions std::vector & fpositions( myvals.getFirstAtomVector() ); - if( fpositions.size()!=ablocks.size() ) fpositions.resize( ablocks.size() ); - for(unsigned i=0; i::performTask( const unsigned& task_index, MultiValue second=first+pbcDistance(first,second); } } - } else if( fpositions.size()==1 ) fpositions[0]=delta(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) ); + } else if( fpositions.size()==1 ) { + fpositions[0]=delta(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) ); + } // Retrieve the masses and charges myvals.resizeTemporyVector(2); std::vector & mass( myvals.getTemporyVector(0) ); std::vector & charge( myvals.getTemporyVector(1) ); - if( mass.size()!=ablocks.size() ) { mass.resize(ablocks.size()); charge.resize(ablocks.size()); } - for(unsigned i=0; i values( getNumberOfComponents() ); std::vector & virial( myvals.getFirstAtomVirialVector() ); std::vector > & derivs( myvals.getFirstAtomDerivativeVector() ); - if( derivs.size()!=values.size() ) { derivs.resize( values.size() ); virial.resize( values.size() ); } + if( derivs.size()!=values.size() ) { + derivs.resize( values.size() ); + virial.resize( values.size() ); + } for(unsigned i=0; i& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); + // active methods: void calculate() override; - static void calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); + MULTICOLVAR_SETTINGS(multiColvars::emptyMode); }; typedef ColvarShortcut PlaneShortcut; @@ -104,11 +102,11 @@ void Plane::parseAtomList( const int& num, std::vector& atoms, Actio } else if( num<0 || atoms.size()>0 ) aa->error("Number of specified atoms should be either 3 or 4"); } -unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) { +Plane::Modetype Plane::getModeAndSetupValues( ActionWithValue* av ) { av->addComponentWithDerivatives("x"); av->componentIsNotPeriodic("x"); av->addComponentWithDerivatives("y"); av->componentIsNotPeriodic("y"); av->addComponentWithDerivatives("z"); av->componentIsNotPeriodic("z"); - return 0; + return {}; } Plane::Plane(const ActionOptions&ao): @@ -127,7 +125,7 @@ Plane::Plane(const ActionOptions&ao): if(pbc) log.printf(" using periodic boundary conditions\n"); else log.printf(" without periodic boundary conditions\n"); - unsigned mode = getModeAndSetupValues( this ); + /*Modetype mode = */getModeAndSetupValues( this ); requestAtoms(atoms); checkRead(); } @@ -135,13 +133,13 @@ Plane::Plane(const ActionOptions&ao): void Plane::calculate() { if(pbc) makeWhole(); - calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this ); + calculateCV( {}, masses, charges, getPositions(), value, derivs, virial, this ); setValue( value[0] ); for(unsigned i=0; i& masses, const std::vector& charges, +void Plane::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) { Vector d1=delta( pos[1], pos[0] ); diff --git a/src/colvar/Position.cpp b/src/colvar/Position.cpp index 20bdc8f234..b9cf37ce10 100644 --- a/src/colvar/Position.cpp +++ b/src/colvar/Position.cpp @@ -100,12 +100,9 @@ class Position : public Colvar { static void registerKeywords( Keywords& keys ); explicit Position(const ActionOptions&); static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); // active methods: void calculate() override; - static void calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); + MULTICOLVAR_SETTINGS(multiColvars::plainOrScaled); }; typedef ColvarShortcut PositionShortcut; @@ -139,8 +136,10 @@ Position::Position(const ActionOptions&ao): { for(unsigned i=0; i<3; ++i) derivs[i].resize(1); std::vector atoms; parseAtomList(-1,atoms,this); - unsigned mode=getModeAndSetupValues(this); - if( mode==1 ) scaled_components=true; + Modetype mode=getModeAndSetupValues(this); + if( mode==Modetype::scaled ) { + scaled_components=true; + } bool nopbc=!pbc; parseFlag("NOPBC",nopbc); @@ -159,19 +158,20 @@ void Position::parseAtomList( const int& num, std::vector& t, Action else if( num<0 || t.size()!=0 ) aa->error("Number of specified atoms should be 1"); } -unsigned Position::getModeAndSetupValues( ActionWithValue* av ) { - bool sc; av->parseFlag("SCALED_COMPONENTS",sc); +Position::Modetype Position::getModeAndSetupValues( ActionWithValue* av ) { + bool sc; + av->parseFlag("SCALED_COMPONENTS",sc); if(sc) { av->addComponentWithDerivatives("a"); av->componentIsPeriodic("a","-0.5","+0.5"); av->addComponentWithDerivatives("b"); av->componentIsPeriodic("b","-0.5","+0.5"); av->addComponentWithDerivatives("c"); av->componentIsPeriodic("c","-0.5","+0.5"); - return 1; + return Modetype::scaled; } av->addComponentWithDerivatives("x"); av->componentIsNotPeriodic("x"); av->addComponentWithDerivatives("y"); av->componentIsNotPeriodic("y"); av->addComponentWithDerivatives("z"); av->componentIsNotPeriodic("z"); av->log<<" WARNING: components will not have the proper periodicity - see manual\n"; - return 0; + return Modetype::plain; } // calculator @@ -185,7 +185,7 @@ void Position::calculate() { } if(scaled_components) { - calculateCV( 1, masses, charges, distance, value, derivs, virial, this ); + calculateCV( Modetype::scaled, masses, charges, distance, value, derivs, virial, this ); Value* valuea=getPntrToComponent("a"); Value* valueb=getPntrToComponent("b"); Value* valuec=getPntrToComponent("c"); @@ -196,7 +196,7 @@ void Position::calculate() { setAtomsDerivatives (valuec,0,derivs[2][0]); valuec->set(value[2]); } else { - calculateCV( 0, masses, charges, distance, value, derivs, virial, this ); + calculateCV( Modetype::plain, masses, charges, distance, value, derivs, virial, this ); Value* valuex=getPntrToComponent("x"); Value* valuey=getPntrToComponent("y"); Value* valuez=getPntrToComponent("z"); @@ -215,10 +215,10 @@ void Position::calculate() { } } -void Position::calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, +void Position::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) { - if( mode==1 ) { + if( mode==Modetype::scaled ) { Vector d=aa->getPbc().realToScaled(pos[0]); vals[0]=Tools::pbc(d[0]); vals[1]=Tools::pbc(d[1]); vals[2]=Tools::pbc(d[2]); derivs[0][0]=matmul(aa->getPbc().getInvBox(),Vector(+1,0,0)); diff --git a/src/colvar/SelectMassCharge.cpp b/src/colvar/SelectMassCharge.cpp index 4d5a4855ea..53039e9f26 100644 --- a/src/colvar/SelectMassCharge.cpp +++ b/src/colvar/SelectMassCharge.cpp @@ -87,12 +87,9 @@ class SelectMassCharge : public Colvar { static void registerKeywords( Keywords& keys ); explicit SelectMassCharge(const ActionOptions&); static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); // active methods: void calculate() override; - static void calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); + MULTICOLVAR_SETTINGS(multiColvars::emptyMode); }; typedef ColvarShortcut MQShortcut; @@ -118,7 +115,7 @@ SelectMassCharge::SelectMassCharge(const ActionOptions&ao): PLUMED_COLVAR_INIT(ao) { std::vector atoms; parseAtomList(-1,atoms,this); - unsigned mode=getModeAndSetupValues(this); + /*Modetype mode=*/getModeAndSetupValues(this); requestAtoms(atoms); } @@ -128,7 +125,7 @@ void SelectMassCharge::parseAtomList( const int& num, std::vector& t else if( num<0 || t.size()!=0 ) aa->error("Number of specified atoms should be 1"); } -unsigned SelectMassCharge::getModeAndSetupValues( ActionWithValue* av ) { +SelectMassCharge::Modetype SelectMassCharge::getModeAndSetupValues( ActionWithValue* av ) { av->addValueWithDerivatives(); av->setNotPeriodic(); bool constant=true; ActionAtomistic* aa=dynamic_cast( av ); plumed_assert( aa ); for(unsigned i=0; igetNumberOfAtoms(); ++i) { @@ -138,21 +135,24 @@ unsigned SelectMassCharge::getModeAndSetupValues( ActionWithValue* av ) { } if( !constant ) av->error("cannot deal with non-constant " + av->getName() + " values"); (av->copyOutput(0))->setConstant(); - return 0; + return {}; } // calculator void SelectMassCharge::calculate() { std::vector masses(1), charges(1), vals(1); std::vector pos; std::vector > derivs; std::vector virial; - calculateCV( 0, masses, charges, pos, vals, derivs, virial, this ); setValue( vals[0] ); + calculateCV( {}, masses, charges, pos, vals, derivs, virial, this ); setValue( vals[0] ); } -void SelectMassCharge::calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, +void SelectMassCharge::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) { - if( aa->getName().find("MASSES")!=std::string::npos ) vals[0]=masses[0]; - else if( aa->chargesWereSet ) vals[0]=charges[0]; + if( aa->getName().find("MASSES")!=std::string::npos ) { + vals[0]=masses[0]; + } else if( aa->chargesWereSet ) { + vals[0]=charges[0]; + } } } diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp index 710b1359eb..e22848390d 100644 --- a/src/colvar/Torsion.cpp +++ b/src/colvar/Torsion.cpp @@ -116,13 +116,13 @@ class Torsion : public Colvar { public: explicit Torsion(const ActionOptions&); static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); // active methods: void calculate() override; - static void calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); static void registerKeywords(Keywords& keys); + enum class torsionModes { + torsion,cosine + }; + MULTICOLVAR_SETTINGS(torsionModes); }; typedef ColvarShortcut TorsionShortcut; @@ -168,8 +168,10 @@ Torsion::Torsion(const ActionOptions&ao): log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n", v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial()); } else parseAtomList(-1,atoms,this); - unsigned mode=getModeAndSetupValues(this); - if( mode==1 ) do_cosine=true; + Modetype mode=getModeAndSetupValues(this); + if( mode==Modetype::cosine ) { + do_cosine=true; + } bool nopbc=!pbc; parseFlag("NOPBC",nopbc); @@ -213,25 +215,37 @@ void Torsion::parseAtomList( const int& num, std::vector& t, ActionA } else if( t.size()!=4 ) aa->error("ATOMS should specify 4 atoms"); } -unsigned Torsion::getModeAndSetupValues( ActionWithValue* av ) { - bool do_cos; av->parseFlag("COSINE",do_cos); - if(do_cos) av->log.printf(" calculating cosine instead of torsion\n"); - +Torsion::Modetype Torsion::getModeAndSetupValues( ActionWithValue* av ) { + bool do_cos; + av->parseFlag("COSINE",do_cos); av->addValueWithDerivatives(); - if(!do_cos) { av->setPeriodic("-pi","pi"); return 0; } - av->setNotPeriodic(); return 1; + if(do_cos) { + av->log.printf(" calculating cosine instead of torsion\n"); + av->setNotPeriodic(); + return Modetype::cosine; + } + av->setPeriodic("-pi","pi"); + return Modetype::torsion; } // calculator void Torsion::calculate() { - if(pbc) makeWhole(); - if(do_cosine) calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this ); - else calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this ); - for(unsigned i=0; i<6; ++i) setAtomsDerivatives(i,derivs[0][i] ); - setValue(value[0]); setBoxDerivatives( virial[0] ); + if(pbc) { + makeWhole(); + } + if(do_cosine) { + calculateCV( Modetype::cosine, masses, charges, getPositions(), value, derivs, virial, this ); + } else { + calculateCV( Modetype::torsion, masses, charges, getPositions(), value, derivs, virial, this ); + } + for(unsigned i=0; i<6; ++i) { + setAtomsDerivatives(i,derivs[0][i] ); + } + setValue(value[0]); + setBoxDerivatives( virial[0] ); } -void Torsion::calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, +void Torsion::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) { Vector d0=delta(pos[1],pos[0]); @@ -240,7 +254,7 @@ void Torsion::calculateCV( const unsigned& mode, const std::vector& mass Vector dd0,dd1,dd2; PLMD::Torsion t; vals[0] = t.compute(d0,d1,d2,dd0,dd1,dd2); - if(mode==1) { + if(mode==Modetype::cosine) { dd0 *= -std::sin(vals[0]); dd1 *= -std::sin(vals[0]); dd2 *= -std::sin(vals[0]); diff --git a/src/core/ActionWithVector.cpp b/src/core/ActionWithVector.cpp index 8eed7e7771..2bc7f7517a 100644 --- a/src/core/ActionWithVector.cpp +++ b/src/core/ActionWithVector.cpp @@ -277,12 +277,16 @@ void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) { void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) const { const ActionWithMatrix* am = dynamic_cast(this); - myvals.setTaskIndex(current); myvals.vector_call=true; performTask( current, myvals ); + myvals.setTaskIndex(current); + myvals.vector_call=true; + performTask( current, myvals ); for(unsigned i=0; ihasDerivatives() ) continue; + if( am || myval->hasDerivatives() ) + continue; Value* myv = const_cast( myval ); - if( getName()=="RMSD_VECTOR" && myv->getRank()==2 ) continue; + if( getName()=="RMSD_VECTOR" && myv->getRank()==2 ) + continue; myv->set( current, myvals.get( i ) ); } } diff --git a/src/core/Value.cpp b/src/core/Value.cpp index 18fe337d15..68b71ac860 100644 --- a/src/core/Value.cpp +++ b/src/core/Value.cpp @@ -361,7 +361,7 @@ void Value::print( OFile& ofile ) const { void Value::printForce( OFile& ofile ) const { if( shape.size()==0 || getNumberOfValues()==1 ) { - ofile.printField( name, getForce(0) ); + ofile.printField( name, getForce(0) ); } else { std::vector indices( shape.size() ); for(unsigned i=0; i& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); + // active methods: void calculate() override; - static void calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ); + MULTICOLVAR_SETTINGS(::PLMD::colvar::multiColvars::emptyMode); }; typedef colvar::ColvarShortcut QuaternionShortcut; @@ -145,7 +143,7 @@ Quaternion::Quaternion(const ActionOptions&ao): parseFlag("NOPBC",nopbc); pbc=!nopbc; - unsigned mode = getModeAndSetupValues( this ); + /*Modetype mode =*/ getModeAndSetupValues( this ); requestAtoms(atoms); } @@ -154,19 +152,19 @@ void Quaternion::parseAtomList( const int& num, std::vector& t, Acti if( t.size()==3 ) aa->log.printf(" involving atoms %d %d %d\n",t[0].serial(),t[1].serial(),t[0].serial()); } -unsigned Quaternion::getModeAndSetupValues( ActionWithValue* av ) { +Quaternion::Modetype Quaternion::getModeAndSetupValues( ActionWithValue* av ) { // This sets up values that we can pass around in PLUMED av->addComponentWithDerivatives("w"); av->componentIsNotPeriodic("w"); av->addComponentWithDerivatives("i"); av->componentIsNotPeriodic("i"); av->addComponentWithDerivatives("j"); av->componentIsNotPeriodic("j"); av->addComponentWithDerivatives("k"); av->componentIsNotPeriodic("k"); - return 0; + return {}; } void Quaternion::calculate() { if(pbc) makeWhole(); - calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this ); + calculateCV( {}, masses, charges, getPositions(), value, derivs, virial, this ); for(unsigned j=0; j<4; ++j) { Value* valuej=getPntrToComponent(j); for(unsigned i=0; i<3; ++i) setAtomsDerivatives(valuej,i,derivs[j][i] ); @@ -176,7 +174,7 @@ void Quaternion::calculate() { } // calculator -void Quaternion::calculateCV( const unsigned& mode, const std::vector& masses, const std::vector& charges, +void Quaternion::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) { //declarations From 7957a2eefd030a1eb4fef301b9006679b267f6d2 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 20 Nov 2024 08:48:50 +0100 Subject: [PATCH 02/10] parseAtomList is a MULTICOLVAR needed function --- src/colvar/Angle.cpp | 5 ++--- src/colvar/DihedralCorrelation.cpp | 5 ++--- src/colvar/Dipole.cpp | 6 ++--- src/colvar/Distance.cpp | 5 ++--- src/colvar/MultiColvarTemplate.h | 36 ++++++++++++++++++++---------- src/colvar/Plane.cpp | 6 ++--- src/colvar/Position.cpp | 5 ++--- src/colvar/SelectMassCharge.cpp | 5 ++--- src/colvar/Torsion.cpp | 5 ++--- src/crystdistrib/Quaternion.cpp | 2 +- 10 files changed, 41 insertions(+), 39 deletions(-) diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp index 0a6579cd00..ac005ac0fa 100644 --- a/src/colvar/Angle.cpp +++ b/src/colvar/Angle.cpp @@ -104,12 +104,11 @@ class Angle : public Colvar { std::vector > derivs; std::vector virial; public: - MULTICOLVAR_SETTINGS(multiColvars::emptyMode); + MULTICOLVAR_DEFAULT(multiColvars::emptyMode); explicit Angle(const ActionOptions&); // active methods: void calculate() override; static void registerKeywords( Keywords& keys ); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); }; @@ -126,7 +125,7 @@ void Angle::registerKeywords( Keywords& keys ) { keys.setValueDescription("the ANGLE involving these atoms"); } -void Angle::parseAtomList( const int& num, std::vector& atoms, ActionAtomistic* aa ) { +void Angle::parseAtomList( int const num, std::vector& atoms, ActionAtomistic* aa ) { aa->parseAtomList("ATOMS",num,atoms); if(atoms.size()==3) { aa->log.printf(" between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial()); diff --git a/src/colvar/DihedralCorrelation.cpp b/src/colvar/DihedralCorrelation.cpp index 28ff787af9..913e661594 100644 --- a/src/colvar/DihedralCorrelation.cpp +++ b/src/colvar/DihedralCorrelation.cpp @@ -68,10 +68,9 @@ class DihedralCorrelation : public Colvar { std::vector > derivs; std::vector virial; public: - MULTICOLVAR_SETTINGS(multiColvars::emptyMode); + MULTICOLVAR_DEFAULT(multiColvars::emptyMode); static void registerKeywords( Keywords& keys ); explicit DihedralCorrelation(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); void calculate() override; }; @@ -107,7 +106,7 @@ DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao): else log.printf(" without periodic boundary conditions\n"); } -void DihedralCorrelation::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { +void DihedralCorrelation::parseAtomList( int const num, std::vector& t, ActionAtomistic* aa ) { aa->parseAtomList("ATOMS",num,t); if( num<0 && t.size()!=8 ) aa->error("Number of specified atoms should be 8"); if( t.size()==8 ) { diff --git a/src/colvar/Dipole.cpp b/src/colvar/Dipole.cpp index 10873ca79c..a42de3f65f 100644 --- a/src/colvar/Dipole.cpp +++ b/src/colvar/Dipole.cpp @@ -90,11 +90,9 @@ class Dipole : public Colvar { Value* valuez=nullptr; public: explicit Dipole(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); void calculate() override; static void registerKeywords(Keywords& keys); - MULTICOLVAR_SETTINGS_MODE(multiColvars::components); - MULTICOLVAR_SETTINGS_SETUPF(); + MULTICOLVAR_SETTINGS_BASE(multiColvars::components); //Declaring this playinly becasue we are using modifying the charges static void calculateCV( Modetype mode, const std::vector& masses, @@ -149,7 +147,7 @@ Dipole::Dipole(const ActionOptions&ao): requestAtoms(ga_lista); } -void Dipole::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { +void Dipole::parseAtomList( int num, std::vector& t, ActionAtomistic* aa ) { aa->parseAtomList("GROUP",num,t); if( t.size()>0 ) { aa->log.printf(" of %u atoms\n",static_cast(t.size())); diff --git a/src/colvar/Distance.cpp b/src/colvar/Distance.cpp index 7145707e36..6f5ad0e551 100644 --- a/src/colvar/Distance.cpp +++ b/src/colvar/Distance.cpp @@ -134,10 +134,9 @@ class Distance : public Colvar { public: static void registerKeywords( Keywords& keys ); explicit Distance(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); // active methods: void calculate() override; - MULTICOLVAR_SETTINGS(multiColvars::scaledComponents); + MULTICOLVAR_DEFAULT(multiColvars::scaledComponents); }; typedef ColvarShortcut DistanceShortcut; @@ -199,7 +198,7 @@ Distance::Distance(const ActionOptions&ao): requestAtoms(atoms); } -void Distance::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { +void Distance::parseAtomList( int const num, std::vector& t, ActionAtomistic* aa ) { aa->parseAtomList("ATOMS",num,t); if( t.size()==2 ) aa->log.printf(" between atoms %d %d\n",t[0].serial(),t[1].serial()); } diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h index 578aedcbc7..998ae5cc77 100644 --- a/src/colvar/MultiColvarTemplate.h +++ b/src/colvar/MultiColvarTemplate.h @@ -43,8 +43,9 @@ enum class scaledComponents { noComponents }; } -#define MULTICOLVAR_SETTINGS_MODE(type) using Modetype=type -#define MULTICOLVAR_SETTINGS_SETUPF() static Modetype getModeAndSetupValues( ActionWithValue* av ) +#define MULTICOLVAR_SETTINGS_BASE(type) using Modetype=type; \ + static void parseAtomList( int num, std::vector& t, ActionAtomistic* aa ); \ + static Modetype getModeAndSetupValues( ActionWithValue* av ) #define MULTICOLVAR_SETTINGS_CALCULATE_CONST() static void calculateCV( Modetype mode, \ const std::vector& masses, \ const std::vector& charges, \ @@ -54,8 +55,7 @@ enum class scaledComponents { std::vector& virial, \ const ActionAtomistic* aa ) -#define MULTICOLVAR_SETTINGS(type) MULTICOLVAR_SETTINGS_MODE(type); \ - MULTICOLVAR_SETTINGS_SETUPF(); \ +#define MULTICOLVAR_DEFAULT(type) MULTICOLVAR_SETTINGS_BASE(type); \ MULTICOLVAR_SETTINGS_CALCULATE_CONST() @@ -106,34 +106,46 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): std::vector all_atoms; if( getName()=="POSITION_VECTOR" || getName()=="MASS_VECTOR" || getName()=="CHARGE_VECTOR" ) parseAtomList( "ATOMS", all_atoms ); if( all_atoms.size()>0 ) { - ablocks.resize(1); ablocks[0].resize( all_atoms.size() ); - for(unsigned i=0; i t; for(int i=1;; ++i ) { Colvar::parseAtomList( i, t, this ); if( t.empty() ) break; - if( i==1 ) { ablocks.resize(t.size()); } + if( i==1 ) { + ablocks.resize(t.size()); + } if( t.size()!=ablocks.size() ) { - std::string ss; Tools::convert(i,ss); + std::string ss; + Tools::convert(i,ss); error("ATOMS" + ss + " keyword has the wrong number of atoms"); } for(unsigned j=0; j& t, ActionAtomistic* aa ); - // active methods: void calculate() override; - MULTICOLVAR_SETTINGS(multiColvars::emptyMode); + MULTICOLVAR_DEFAULT(multiColvars::emptyMode); }; typedef ColvarShortcut PlaneShortcut; @@ -90,7 +88,7 @@ void Plane::registerKeywords( Keywords& keys ) { keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); } -void Plane::parseAtomList( const int& num, std::vector& atoms, ActionAtomistic* aa ) { +void Plane::parseAtomList( int const num, std::vector& atoms, ActionAtomistic* aa ) { aa->parseAtomList("ATOMS",num,atoms); if(atoms.size()==3) { aa->log.printf(" containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial()); diff --git a/src/colvar/Position.cpp b/src/colvar/Position.cpp index b9cf37ce10..53dd180bbb 100644 --- a/src/colvar/Position.cpp +++ b/src/colvar/Position.cpp @@ -99,10 +99,9 @@ class Position : public Colvar { public: static void registerKeywords( Keywords& keys ); explicit Position(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); // active methods: void calculate() override; - MULTICOLVAR_SETTINGS(multiColvars::plainOrScaled); + MULTICOLVAR_DEFAULT(multiColvars::plainOrScaled); }; typedef ColvarShortcut PositionShortcut; @@ -152,7 +151,7 @@ Position::Position(const ActionOptions&ao): requestAtoms(atoms); } -void Position::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { +void Position::parseAtomList(int const num, std::vector& t, ActionAtomistic* aa ) { aa->parseAtomList("ATOM",num,t); if( t.size()==1 ) aa->log.printf(" for atom %d\n",t[0].serial()); else if( num<0 || t.size()!=0 ) aa->error("Number of specified atoms should be 1"); diff --git a/src/colvar/SelectMassCharge.cpp b/src/colvar/SelectMassCharge.cpp index 53039e9f26..e3e98cccac 100644 --- a/src/colvar/SelectMassCharge.cpp +++ b/src/colvar/SelectMassCharge.cpp @@ -86,10 +86,9 @@ class SelectMassCharge : public Colvar { public: static void registerKeywords( Keywords& keys ); explicit SelectMassCharge(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); // active methods: void calculate() override; - MULTICOLVAR_SETTINGS(multiColvars::emptyMode); + MULTICOLVAR_DEFAULT(multiColvars::emptyMode); }; typedef ColvarShortcut MQShortcut; @@ -119,7 +118,7 @@ SelectMassCharge::SelectMassCharge(const ActionOptions&ao): requestAtoms(atoms); } -void SelectMassCharge::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { +void SelectMassCharge::parseAtomList( int const num, std::vector& t, ActionAtomistic* aa ) { aa->parseAtomList("ATOM",num,t); if( t.size()==1 ) aa->log.printf(" for atom %d\n",t[0].serial()); else if( num<0 || t.size()!=0 ) aa->error("Number of specified atoms should be 1"); diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp index e22848390d..2efc5f2b2e 100644 --- a/src/colvar/Torsion.cpp +++ b/src/colvar/Torsion.cpp @@ -115,14 +115,13 @@ class Torsion : public Colvar { std::vector virial; public: explicit Torsion(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); // active methods: void calculate() override; static void registerKeywords(Keywords& keys); enum class torsionModes { torsion,cosine }; - MULTICOLVAR_SETTINGS(torsionModes); + MULTICOLVAR_DEFAULT(torsionModes); }; typedef ColvarShortcut TorsionShortcut; @@ -183,7 +182,7 @@ Torsion::Torsion(const ActionOptions&ao): requestAtoms(atoms); } -void Torsion::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { +void Torsion::parseAtomList( int const num, std::vector& t, ActionAtomistic* aa ) { std::vector v1,v2,axis; aa->parseAtomList("ATOMS",num,t); aa->parseAtomList("VECTORA",num,v1); diff --git a/src/crystdistrib/Quaternion.cpp b/src/crystdistrib/Quaternion.cpp index 574ef3524c..4029abd40d 100644 --- a/src/crystdistrib/Quaternion.cpp +++ b/src/crystdistrib/Quaternion.cpp @@ -107,7 +107,7 @@ class Quaternion : public Colvar { // active methods: void calculate() override; - MULTICOLVAR_SETTINGS(::PLMD::colvar::multiColvars::emptyMode); + MULTICOLVAR_DEFAULT(::PLMD::colvar::multiColvars::emptyMode); }; typedef colvar::ColvarShortcut QuaternionShortcut; From afbbf8229335c50d12133235149db957d4a610ad Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 20 Nov 2024 09:21:28 +0100 Subject: [PATCH 03/10] Writing a guide for myself --- src/colvar/MultiColvarTemplate.h | 30 +++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h index 998ae5cc77..7f10b4ec36 100644 --- a/src/colvar/MultiColvarTemplate.h +++ b/src/colvar/MultiColvarTemplate.h @@ -43,6 +43,30 @@ enum class scaledComponents { noComponents }; } +/*MANUAL DRAFT +To setup a CV compatible with multicolvartemplate you need to add + - A type that will be used to pass to calculateCV the calculation settings + - this type will be called `Modetype` + - using Modetype=type; + - some default types are already defined in MultiColvarTemplate.h + - for example `using Modetype=PLME::colvar::multiColvars::emptyMode;` when no setting inmact the calculation + - others predefined types are: + - `enum class components {withCompontents, noComponents};` + - `enum class plainOrScaled {scaled, plain};` + - `enum class scaledComponents {withCompontents, scaledComponents, noComponents};` + - some static function that will be called by the MultiColvarTemplate: + - these functions will need to match the foloowing signatures (withour the constness of the non const& arguments): + - `static void parseAtomList( int num, std::vector& t, ActionAtomistic* aa );` + - `Modetype getModeAndSetupValues( ActionWithValue* av )` + - `static void calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) + - this function will be called by MultiColvarTemplate to calculate the CV value on the inputs + - to avoid code repetitition you should change ::calculate() to call this function + - By default all the inputs are const ref, but the constedness can be changed, since the MulticolvarTemplate will pass the plain references + - By default you can decorate the public part of your CV wiht `MULTICOLVAR_DEFAULT(modetype here)` to "conver" the declaration of the CV in MultiColvarTemplate-compatible one, you will need to manually declare and implement the methods + - If you need to change some constness of the `calculateCV` signature use `MULTICOLVAR_SETTINGS_BASE(modetype here)` and then declare `calculateCV` with the desired constness in the signature (see Dipole.cpp) +*/ + + #define MULTICOLVAR_SETTINGS_BASE(type) using Modetype=type; \ static void parseAtomList( int num, std::vector& t, ActionAtomistic* aa ); \ static Modetype getModeAndSetupValues( ActionWithValue* av ) @@ -119,7 +143,7 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): if( i==1 ) { ablocks.resize(t.size()); - } + } if( t.size()!=ablocks.size() ) { std::string ss; Tools::convert(i,ss); @@ -132,7 +156,7 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): t.resize(0); } } - if( all_atoms.size()==0 ) { + if( all_atoms.size()==0 ){ error("No atoms have been specified"); } requestAtoms(all_atoms); @@ -143,7 +167,7 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): } if( keywords.exists("WHOLEMOLECULES") ) { parseFlag("WHOLEMOLECULES",wholemolecules); - if( wholemolecules ) { + if( wholemolecules ){ usepbc=false; } } From 8aee146412b06beae6ccab3f2a840b9684a0ccaf Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 20 Nov 2024 10:52:50 +0100 Subject: [PATCH 04/10] Using swithches in distance, less copy-paste :) --- src/colvar/Distance.cpp | 97 ++++++++++++++++++-------------- src/colvar/MultiColvarTemplate.h | 6 +- 2 files changed, 59 insertions(+), 44 deletions(-) diff --git a/src/colvar/Distance.cpp b/src/colvar/Distance.cpp index 6f5ad0e551..05f3246c55 100644 --- a/src/colvar/Distance.cpp +++ b/src/colvar/Distance.cpp @@ -124,19 +124,20 @@ Calculate a vector containing the distances between various pairs of atoms //+ENDPLUMEDOC class Distance : public Colvar { - bool components; - bool scaled_components; - bool pbc; - - std::vector value, masses, charges; - std::vector > derivs; - std::vector virial; public: static void registerKeywords( Keywords& keys ); explicit Distance(const ActionOptions&); // active methods: void calculate() override; MULTICOLVAR_DEFAULT(multiColvars::scaledComponents); +private: + std::vector value{{0.0}}; + std::vectormasses{{0.0}}; + std::vector charges{{0.0}}; + std::vector > derivs{std::vector{Vector{},Vector{}}}; + std::vector virial{Tensor()}; + Modetype mode_; + bool pbc=true; }; typedef ColvarShortcut DistanceShortcut; @@ -161,13 +162,7 @@ void Distance::registerKeywords( Keywords& keys ) { } Distance::Distance(const ActionOptions&ao): - PLUMED_COLVAR_INIT(ao), - components(false), - scaled_components(false), - pbc(true), - value(1), - derivs(1), - virial(1) + PLUMED_COLVAR_INIT(ao) { derivs[0].resize(2); std::vector atoms; @@ -182,13 +177,8 @@ Distance::Distance(const ActionOptions&ao): if(pbc) log.printf(" using periodic boundary conditions\n"); else log.printf(" without periodic boundary conditions\n"); - Modetype mode = getModeAndSetupValues( this ); - if(mode==Modetype::withCompontents) { - components=true; - } else if(mode==Modetype::scaledComponents) { - scaled_components=true; - } - if( components || scaled_components ) { + mode_ = getModeAndSetupValues( this ); + if(mode_==Modetype::withCompontents || mode_==Modetype::scaledComponents) { value.resize(3); derivs.resize(3); virial.resize(3); @@ -231,42 +221,61 @@ Distance::Modetype Distance::getModeAndSetupValues( ActionWithValue* av ) { void Distance::calculate() { if(pbc) makeWhole(); + calculateCV( mode_, masses, charges, getPositions(), value, derivs, virial, this ); - if( components ) { - calculateCV( Modetype::withCompontents, masses, charges, getPositions(), value, derivs, virial, this ); + switch (mode_) { + case Modetype::withCompontents: { Value* valuex=getPntrToComponent("x"); Value* valuey=getPntrToComponent("y"); Value* valuez=getPntrToComponent("z"); - for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuex,i,derivs[0][i] ); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuex,i,derivs[0][i] ); + } setBoxDerivatives(valuex,virial[0]); valuex->set(value[0]); - for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuey,i,derivs[1][i] ); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuey,i,derivs[1][i] ); + } setBoxDerivatives(valuey,virial[1]); valuey->set(value[1]); - for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuez,i,derivs[2][i] ); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuez,i,derivs[2][i] ); + } setBoxDerivatives(valuez,virial[2]); valuez->set(value[2]); - } else if( scaled_components ) { - calculateCV( Modetype::scaledComponents, masses, charges, getPositions(), value, derivs, virial, this ); + } + break; + case Modetype::scaledComponents: { Value* valuea=getPntrToComponent("a"); Value* valueb=getPntrToComponent("b"); Value* valuec=getPntrToComponent("c"); - for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuea,i,derivs[0][i] ); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuea,i,derivs[0][i] ); + } valuea->set(value[0]); - for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valueb,i,derivs[1][i] ); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valueb,i,derivs[1][i] ); + } valueb->set(value[1]); - for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuec,i,derivs[2][i] ); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuec,i,derivs[2][i] ); + } valuec->set(value[2]); - } else { - calculateCV( Modetype::noComponents, masses, charges, getPositions(), value, derivs, virial, this ); - for(unsigned i=0; i<2; ++i) setAtomsDerivatives(i,derivs[0][i] ); + } + break; + + case Modetype::noComponents: { + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(i,derivs[0][i] ); + } setBoxDerivatives(virial[0]); setValue (value[0]); } + } } void Distance::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, @@ -276,7 +285,8 @@ void Distance::calculateCV( Modetype mode, const std::vector& masses, co const double value=distance.modulo(); const double invvalue=1.0/value; - if(mode==Modetype::withCompontents) { + switch (mode) { + case Modetype::withCompontents: { derivs[0][0] = Vector(-1,0,0); derivs[0][1] = Vector(+1,0,0); vals[0] = distance[0]; @@ -289,7 +299,10 @@ void Distance::calculateCV( Modetype mode, const std::vector& masses, co derivs[2][1] = Vector(0,0,+1); vals[2] = distance[2]; setBoxDerivativesNoPbc( pos, derivs, virial ); - } else if(mode==Modetype::scaledComponents) { + } + break; + + case Modetype::scaledComponents: { 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)); @@ -300,16 +313,18 @@ void Distance::calculateCV( Modetype mode, const std::vector& masses, co derivs[2][0] = matmul(aa->getPbc().getInvBox(),Vector(0,0,-1)); derivs[2][1] = matmul(aa->getPbc().getInvBox(),Vector(0,0,+1)); vals[2] = Tools::pbc(d[2]); - } else { + } + break; + + case Modetype::noComponents: { derivs[0][0] = -invvalue*distance; derivs[0][1] = invvalue*distance; setBoxDerivativesNoPbc( pos, derivs, virial ); vals[0] = value; } -} + } } -} - - +} // namespace colvar +} // namespace PLMD diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h index 7f10b4ec36..dacf1e8412 100644 --- a/src/colvar/MultiColvarTemplate.h +++ b/src/colvar/MultiColvarTemplate.h @@ -143,7 +143,7 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): if( i==1 ) { ablocks.resize(t.size()); - } + } if( t.size()!=ablocks.size() ) { std::string ss; Tools::convert(i,ss); @@ -156,7 +156,7 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): t.resize(0); } } - if( all_atoms.size()==0 ){ + if( all_atoms.size()==0 ) { error("No atoms have been specified"); } requestAtoms(all_atoms); @@ -167,7 +167,7 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): } if( keywords.exists("WHOLEMOLECULES") ) { parseFlag("WHOLEMOLECULES",wholemolecules); - if( wholemolecules ){ + if( wholemolecules ) { usepbc=false; } } From 740f6b6466714614b908371c23ede89443cc5a95 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 20 Nov 2024 12:18:29 +0100 Subject: [PATCH 05/10] Simplyfying the signature of calculateCV --- src/colvar/Angle.cpp | 9 ++-- src/colvar/DihedralCorrelation.cpp | 9 ++-- src/colvar/Dipole.cpp | 68 +++++++++++++---------- src/colvar/Distance.cpp | 11 ++-- src/colvar/MultiColvarTemplate.h | 87 +++++++++++++++++++++--------- src/colvar/Plane.cpp | 8 +-- src/colvar/Position.cpp | 56 ++++++++++--------- src/colvar/SelectMassCharge.cpp | 8 +-- src/colvar/Torsion.cpp | 32 +++++------ src/crystdistrib/Quaternion.cpp | 9 ++-- 10 files changed, 181 insertions(+), 116 deletions(-) diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp index ac005ac0fa..f5b9e4dd38 100644 --- a/src/colvar/Angle.cpp +++ b/src/colvar/Angle.cpp @@ -168,7 +168,7 @@ Angle::Angle(const ActionOptions&ao): void Angle::calculate() { if(pbc) makeWhole(); - calculateCV( {}, masses, charges, getPositions(), value, derivs, virial, this ); + calculateCV( {}, masses, charges, getPositions(), multiColvars::Ouput(value, derivs, virial), this ); setValue( value[0] ); for(unsigned i=0; i& /*masses*/, const std::vector& /*charges*/, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* /*aa*/ ) { + const std::vector& pos, + multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); + auto & derivs=out.derivs(); + auto & virial=out.virial(); Vector dij,dik; dij=delta(pos[2],pos[3]); dik=delta(pos[1],pos[0]); diff --git a/src/colvar/DihedralCorrelation.cpp b/src/colvar/DihedralCorrelation.cpp index 913e661594..1a5e42c604 100644 --- a/src/colvar/DihedralCorrelation.cpp +++ b/src/colvar/DihedralCorrelation.cpp @@ -124,15 +124,18 @@ DihedralCorrelation::Modetype DihedralCorrelation::getModeAndSetupValues( Action void DihedralCorrelation::calculate() { if(pbc) makeWhole(); - calculateCV( {}, masses, charges, getPositions(), value, derivs, virial, this ); + calculateCV( {}, masses, charges, getPositions(), multiColvars::Ouput(value, derivs, virial), this ); setValue( value[0] ); for(unsigned i=0; i& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { + const std::vector& pos, + multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); + auto & derivs=out.derivs(); + auto & virial=out.virial(); const Vector d10=delta(pos[1],pos[0]); const Vector d11=delta(pos[2],pos[1]); const Vector d12=delta(pos[3],pos[2]); diff --git a/src/colvar/Dipole.cpp b/src/colvar/Dipole.cpp index a42de3f65f..972fc92b33 100644 --- a/src/colvar/Dipole.cpp +++ b/src/colvar/Dipole.cpp @@ -79,15 +79,6 @@ Calculate a vector of dipole moments for a set of groups of atoms. //+ENDPLUMEDOC class Dipole : public Colvar { - std::vector ga_lista; - bool components; - bool nopbc; - std::vector value, masses, charges; - std::vector > derivs; - std::vector virial; - Value* valuex=nullptr; - Value* valuey=nullptr; - Value* valuez=nullptr; public: explicit Dipole(const ActionOptions&); void calculate() override; @@ -98,10 +89,18 @@ class Dipole : public Colvar { const std::vector& masses, std::vector& charges, const std::vector& pos, - std::vector& vals, - std::vector >& derivs, - std::vector& virial, + multiColvars::Ouput out, const ActionAtomistic* aa ); +private: + std::vector ga_lista; + std::vector value, masses, charges; + std::vector > derivs; + std::vector virial; + Value* valuex=nullptr; + Value* valuey=nullptr; + Value* valuez=nullptr; + Modetype components; + bool nopbc=false; }; typedef ColvarShortcut DipoleShortcut; @@ -123,21 +122,24 @@ void Dipole::registerKeywords(Keywords& keys) { Dipole::Dipole(const ActionOptions&ao): PLUMED_COLVAR_INIT(ao), - components(false), + value(1), derivs(1), - virial(1) + virial(1), + components(getModeAndSetupValues(this)) { parseAtomList(-1,ga_lista,this); charges.resize(ga_lista.size()); - components=(getModeAndSetupValues(this)==Modetype::withCompontents); - if( components ) { + // components=getModeAndSetupValues(this); + if( components == Modetype::withCompontents ) { value.resize(3); derivs.resize(3); virial.resize(3); valuex=getPntrToComponent("x"); valuey=getPntrToComponent("y"); valuez=getPntrToComponent("z"); } - for(unsigned i=0; i& masses, std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { - unsigned N=pos.size(); double ctot=0.; - for(unsigned i=0; i& pos,multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); + auto & derivs=out.derivs(); + auto & virial=out.virial(); + unsigned N=pos.size(); + double ctot=0.; + for(unsigned i=0; i& masses, std: if( mode==Modetype::withCompontents ) { for(unsigned i=0; i& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { + const std::vector& pos, + multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); + auto & derivs=out.derivs(); + auto & virial=out.virial(); Vector distance=delta(pos[0],pos[1]); const double value=distance.modulo(); const double invvalue=1.0/value; - + //the compiler will alert me if I forgot a case :) switch (mode) { case Modetype::withCompontents: { derivs[0][0] = Vector(-1,0,0); diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h index dacf1e8412..fb5d6c386c 100644 --- a/src/colvar/MultiColvarTemplate.h +++ b/src/colvar/MultiColvarTemplate.h @@ -42,7 +42,46 @@ enum class scaledComponents { scaledComponents, noComponents }; -} + + +//TODO: test this +class Ouput +{ + //these are const pointers to non const data,meaning YOU can't change the pointer + std::vector * vals_=nullptr; + std::vector >* derivs_=nullptr; + std::vector * virial_=nullptr; + Ouput()=default; +public: +//const because the data in the class is not changing: we are modifying data pointed by the const pointers + std::vector&vals() const {return *vals_;} + std::vector >&derivs() const {return *derivs_;} + std::vector&virial() const {return *virial_;} + Ouput(std::vector& vals, + std::vector >& derivs, + std::vector& virial) + : vals_(&vals), derivs_(&derivs), virial_(&virial) {} + Ouput(Ouput const& other) : vals_(other.vals_), derivs_(other.derivs_), virial_(other.virial_) {}; + Ouput(Ouput && other) noexcept : Ouput() { + swap(*this, other); + }; + // the input is initialized as copy or move ctor + Ouput& operator=(Ouput other) =delete; + // { + // //Self assignment protection is free[cit] + // swap(*this, other); + // return *this; + // }; + friend void swap(Ouput& a, Ouput& b) noexcept { + std::swap(a.vals_, b.vals_); + std::swap(a.derivs_, b.derivs_); + std::swap(a.virial_, b.virial_); + } + ~Ouput() {}//this does not own anything +}; + +} // namespace multiColvars + /*MANUAL DRAFT To setup a CV compatible with multicolvartemplate you need to add - A type that will be used to pass to calculateCV the calculation settings @@ -58,7 +97,7 @@ To setup a CV compatible with multicolvartemplate you need to add - these functions will need to match the foloowing signatures (withour the constness of the non const& arguments): - `static void parseAtomList( int num, std::vector& t, ActionAtomistic* aa );` - `Modetype getModeAndSetupValues( ActionWithValue* av )` - - `static void calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, const std::vector& pos, std::vector& vals, std::vector >& derivs, std::vector& virial, const ActionAtomistic* aa ) + - `static void calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, const std::vector& pos, PLMD::colvars::multiColvars::Ouput out, const ActionAtomistic* aa ) - this function will be called by MultiColvarTemplate to calculate the CV value on the inputs - to avoid code repetitition you should change ::calculate() to call this function - By default all the inputs are const ref, but the constedness can be changed, since the MulticolvarTemplate will pass the plain references @@ -74,9 +113,7 @@ To setup a CV compatible with multicolvartemplate you need to add const std::vector& masses, \ const std::vector& charges, \ const std::vector& pos, \ - std::vector& vals, \ - std::vector >& derivs,\ - std::vector& virial, \ + PLMD::colvar::multiColvars::Ouput out, \ const ActionAtomistic* aa ) #define MULTICOLVAR_DEFAULT(type) MULTICOLVAR_SETTINGS_BASE(type); \ @@ -86,10 +123,10 @@ To setup a CV compatible with multicolvartemplate you need to add //^no ';' here so that the macro will "consume" it and not generate a double semicolon error -template +template class MultiColvarTemplate : public ActionWithVector { private: - using Modetype =typename Colvar::Modetype; + using Modetype =typename CV::Modetype; /// An index that decides what we are calculating Modetype mode; /// Are we using pbc to calculate the CVs @@ -108,9 +145,9 @@ class MultiColvarTemplate : public ActionWithVector { void calculate() override; }; -template -void MultiColvarTemplate::registerKeywords(Keywords& keys ) { - Colvar::registerKeywords( keys ); +template +void MultiColvarTemplate::registerKeywords(Keywords& keys ) { + CV::registerKeywords( keys ); keys.add("optional","MASK","the label for a sparse matrix that should be used to determine which elements of the matrix should be computed"); unsigned nkeys = keys.size(); for(unsigned i=0; i::registerKeywords(Keywords& keys ) { if( keys.outputComponentExists(".#!value") ) keys.setValueDescription("the " + keys.getDisplayName() + " for each set of specified atoms"); } -template -MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): +template +MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): Action(ao), ActionWithVector(ao), mode(), @@ -138,7 +175,7 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): } else { std::vector t; for(int i=1;; ++i ) { - Colvar::parseAtomList( i, t, this ); + CV::parseAtomList( i, t, this ); if( t.empty() ) break; if( i==1 ) { @@ -175,32 +212,32 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): else log.printf(" without periodic boundary conditions\n"); // Setup the values - mode = Colvar::getModeAndSetupValues( this ); + mode = CV::getModeAndSetupValues( this ); } -template -unsigned MultiColvarTemplate::getNumberOfDerivatives() { +template +unsigned MultiColvarTemplate::getNumberOfDerivatives() { return 3*getNumberOfAtoms()+9; } -template -void MultiColvarTemplate::calculate() { +template +void MultiColvarTemplate::calculate() { if( wholemolecules ) makeWhole(); runAllTasks(); } -template -void MultiColvarTemplate::addValueWithDerivatives( const std::vector& shape ) { +template +void MultiColvarTemplate::addValueWithDerivatives( const std::vector& shape ) { std::vector s(1); s[0]=ablocks[0].size(); addValue( s ); } -template -void MultiColvarTemplate::addComponentWithDerivatives( const std::string& name, const std::vector& shape ) { +template +void MultiColvarTemplate::addComponentWithDerivatives( const std::string& name, const std::vector& shape ) { std::vector s(1); s[0]=ablocks[0].size(); addComponent( name, s ); } -template -void MultiColvarTemplate::performTask( const unsigned& task_index, MultiValue& myvals ) const { +template +void MultiColvarTemplate::performTask( const unsigned& task_index, MultiValue& myvals ) const { // Retrieve the positions std::vector & fpositions( myvals.getFirstAtomVector() ); if( fpositions.size()!=ablocks.size() ) { @@ -248,7 +285,7 @@ void MultiColvarTemplate::performTask( const unsigned& task_index, Multi } } // Calculate the CVs using the method in the Colvar - Colvar::calculateCV( mode, mass, charge, fpositions, values, derivs, virial, this ); + CV::calculateCV( mode, mass, charge, fpositions, multiColvars::Ouput{values, derivs, virial}, this ); for(unsigned i=0; i& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { + const std::vector& pos,multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); + auto & derivs=out.derivs(); + auto & virial=out.virial(); Vector d1=delta( pos[1], pos[0] ); Vector d2=delta( pos[2], pos[3] ); Vector cp = crossProduct( d1, d2 ); diff --git a/src/colvar/Position.cpp b/src/colvar/Position.cpp index 53dd180bbb..9e902befd1 100644 --- a/src/colvar/Position.cpp +++ b/src/colvar/Position.cpp @@ -91,17 +91,20 @@ Create a vector that holds the components of the position of a set of atoms. //+ENDPLUMEDOC class Position : public Colvar { - bool scaled_components; - bool pbc; - std::vector value, masses, charges; - std::vector > derivs; - std::vector virial; public: static void registerKeywords( Keywords& keys ); explicit Position(const ActionOptions&); // active methods: void calculate() override; MULTICOLVAR_DEFAULT(multiColvars::plainOrScaled); +private: + Modetype components; + bool pbc; + std::vector value; + std::vector masses; + std::vector charges; + std::vector > derivs; + std::vector virial; }; typedef ColvarShortcut PositionShortcut; @@ -127,7 +130,6 @@ void Position::registerKeywords( Keywords& keys ) { Position::Position(const ActionOptions&ao): PLUMED_COLVAR_INIT(ao), - scaled_components(false), pbc(true), value(3), derivs(3), @@ -135,10 +137,7 @@ Position::Position(const ActionOptions&ao): { for(unsigned i=0; i<3; ++i) derivs[i].resize(1); std::vector atoms; parseAtomList(-1,atoms,this); - Modetype mode=getModeAndSetupValues(this); - if( mode==Modetype::scaled ) { - scaled_components=true; - } + components=getModeAndSetupValues(this); bool nopbc=!pbc; parseFlag("NOPBC",nopbc); @@ -182,9 +181,9 @@ void Position::calculate() { } else { distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0)); } - - if(scaled_components) { - calculateCV( Modetype::scaled, masses, charges, distance, value, derivs, virial, this ); + calculateCV( components, masses, charges, distance, multiColvars::Ouput(value, derivs, virial), this ); + switch (components) { + case Modetype::scaled: { Value* valuea=getPntrToComponent("a"); Value* valueb=getPntrToComponent("b"); Value* valuec=getPntrToComponent("c"); @@ -194,8 +193,10 @@ void Position::calculate() { valueb->set(value[1]); setAtomsDerivatives (valuec,0,derivs[2][0]); valuec->set(value[2]); - } else { - calculateCV( Modetype::plain, masses, charges, distance, value, derivs, virial, this ); + } + break; + + case Modetype::plain: { Value* valuex=getPntrToComponent("x"); Value* valuey=getPntrToComponent("y"); Value* valuez=getPntrToComponent("z"); @@ -212,26 +213,33 @@ void Position::calculate() { setBoxDerivatives (valuez,virial[2]); valuez->set(value[2]); } + } } void Position::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { - if( mode==Modetype::scaled ) { + const std::vector& pos, + multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); + auto & derivs=out.derivs(); + auto & virial=out.virial(); + switch (mode) + { + case Modetype::scaled: { Vector d=aa->getPbc().realToScaled(pos[0]); vals[0]=Tools::pbc(d[0]); vals[1]=Tools::pbc(d[1]); vals[2]=Tools::pbc(d[2]); derivs[0][0]=matmul(aa->getPbc().getInvBox(),Vector(+1,0,0)); derivs[1][0]=matmul(aa->getPbc().getInvBox(),Vector(0,+1,0)); derivs[2][0]=matmul(aa->getPbc().getInvBox(),Vector(0,0,+1)); - } else { + } + break; + + case Modetype::plain: { for(unsigned i=0; i<3; ++i) vals[i]=pos[0][i]; derivs[0][0]=Vector(+1,0,0); derivs[1][0]=Vector(0,+1,0); derivs[2][0]=Vector(0,0,+1); virial[0]=Tensor(pos[0],Vector(-1,0,0)); virial[1]=Tensor(pos[0],Vector(0,-1,0)); virial[2]=Tensor(pos[0],Vector(0,0,-1)); } + } } -} -} - - - +} // namespace colvar +} // namespacs PLMD diff --git a/src/colvar/SelectMassCharge.cpp b/src/colvar/SelectMassCharge.cpp index e3e98cccac..82f417df14 100644 --- a/src/colvar/SelectMassCharge.cpp +++ b/src/colvar/SelectMassCharge.cpp @@ -139,14 +139,14 @@ SelectMassCharge::Modetype SelectMassCharge::getModeAndSetupValues( ActionWithVa // calculator void SelectMassCharge::calculate() { - std::vector masses(1), charges(1), vals(1); + std::vector masses(1), charges(1), value(1); std::vector pos; std::vector > derivs; std::vector virial; - calculateCV( {}, masses, charges, pos, vals, derivs, virial, this ); setValue( vals[0] ); + calculateCV( {}, masses, charges, pos, multiColvars::Ouput(value, derivs, virial), this ); setValue( value[0] ); } void SelectMassCharge::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { + const std::vector& pos,multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); if( aa->getName().find("MASSES")!=std::string::npos ) { vals[0]=masses[0]; } else if( aa->chargesWereSet ) { diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp index 2efc5f2b2e..2a9c6fb314 100644 --- a/src/colvar/Torsion.cpp +++ b/src/colvar/Torsion.cpp @@ -107,12 +107,6 @@ Calculate multiple torsional angles. //+ENDPLUMEDOC class Torsion : public Colvar { - bool pbc; - bool do_cosine; - - std::vector value, masses, charges; - std::vector > derivs; - std::vector virial; public: explicit Torsion(const ActionOptions&); // active methods: @@ -122,6 +116,12 @@ class Torsion : public Colvar { torsion,cosine }; MULTICOLVAR_DEFAULT(torsionModes); +private: + std::vector value, masses, charges; + std::vector > derivs; + std::vector virial; + bool pbc; + Modetype mode; }; typedef ColvarShortcut TorsionShortcut; @@ -146,7 +146,6 @@ void Torsion::registerKeywords(Keywords& keys) { Torsion::Torsion(const ActionOptions&ao): PLUMED_COLVAR_INIT(ao), pbc(true), - do_cosine(false), value(1), derivs(1), virial(1) @@ -167,10 +166,7 @@ Torsion::Torsion(const ActionOptions&ao): log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n", v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial()); } else parseAtomList(-1,atoms,this); - Modetype mode=getModeAndSetupValues(this); - if( mode==Modetype::cosine ) { - do_cosine=true; - } + mode=getModeAndSetupValues(this); bool nopbc=!pbc; parseFlag("NOPBC",nopbc); @@ -232,11 +228,7 @@ void Torsion::calculate() { if(pbc) { makeWhole(); } - if(do_cosine) { - calculateCV( Modetype::cosine, masses, charges, getPositions(), value, derivs, virial, this ); - } else { - calculateCV( Modetype::torsion, masses, charges, getPositions(), value, derivs, virial, this ); - } + calculateCV( mode, masses, charges, getPositions(), multiColvars::Ouput(value, derivs, virial), this ); for(unsigned i=0; i<6; ++i) { setAtomsDerivatives(i,derivs[0][i] ); } @@ -245,8 +237,12 @@ void Torsion::calculate() { } void Torsion::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { + const std::vector& pos, + multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); + auto & derivs=out.derivs(); + auto & virial=out.virial(); + Vector d0=delta(pos[1],pos[0]); Vector d1=delta(pos[3],pos[2]); Vector d2=delta(pos[5],pos[4]); diff --git a/src/crystdistrib/Quaternion.cpp b/src/crystdistrib/Quaternion.cpp index 4029abd40d..86fe8ac21a 100644 --- a/src/crystdistrib/Quaternion.cpp +++ b/src/crystdistrib/Quaternion.cpp @@ -164,7 +164,7 @@ Quaternion::Modetype Quaternion::getModeAndSetupValues( ActionWithValue* av ) { void Quaternion::calculate() { if(pbc) makeWhole(); - calculateCV( {}, masses, charges, getPositions(), value, derivs, virial, this ); + calculateCV( {}, masses, charges, getPositions(), PLMD::colvar::multiColvars::Ouput(value, derivs, virial), this ); for(unsigned j=0; j<4; ++j) { Value* valuej=getPntrToComponent(j); for(unsigned i=0; i<3; ++i) setAtomsDerivatives(valuej,i,derivs[j][i] ); @@ -175,8 +175,11 @@ void Quaternion::calculate() { // calculator void Quaternion::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, - const std::vector& pos, std::vector& vals, std::vector >& derivs, - std::vector& virial, const ActionAtomistic* aa ) { + const std::vector& pos, + PLMD::colvar::multiColvars::Ouput out, const ActionAtomistic* aa ) { + auto & vals=out.vals(); + auto & derivs=out.derivs(); + auto & virial=out.virial(); //declarations Vector vec1_comp = delta( pos[0], pos[1] ); //components between atom 1 and 2 Vector vec2_comp = delta( pos[0], pos[2] ); //components between atom 1 and 3 From 974c3ce99e7d0c61d49ce9153fc7e06146bdb5f9 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 20 Nov 2024 12:56:20 +0100 Subject: [PATCH 06/10] removing the forward declaration of selectMassCharge from ActionAtomistic + scalar mass/charge should work --- src/colvar/SelectMassCharge.cpp | 118 +++++++++++++++++++++++--------- src/core/ActionAtomistic.h | 22 ++++-- 2 files changed, 101 insertions(+), 39 deletions(-) diff --git a/src/colvar/SelectMassCharge.cpp b/src/colvar/SelectMassCharge.cpp index 82f417df14..e396a4c646 100644 --- a/src/colvar/SelectMassCharge.cpp +++ b/src/colvar/SelectMassCharge.cpp @@ -82,6 +82,12 @@ Get the mass of one or multiple atoms namespace PLMD { namespace colvar { +namespace +{ +enum class MC {Mass,Charge}; +} // namespace unnamed + +template class SelectMassCharge : public Colvar { public: static void registerKeywords( Keywords& keys ); @@ -89,18 +95,27 @@ class SelectMassCharge : public Colvar { // active methods: void calculate() override; MULTICOLVAR_DEFAULT(multiColvars::emptyMode); +private: + AtomNumber theAtom; }; -typedef ColvarShortcut MQShortcut; -PLUMED_REGISTER_ACTION(MQShortcut,"MASS") -PLUMED_REGISTER_ACTION(MQShortcut,"CHARGE") -PLUMED_REGISTER_ACTION(SelectMassCharge,"MASS_SCALAR") -PLUMED_REGISTER_ACTION(SelectMassCharge,"CHARGE_SCALAR") -typedef MultiColvarTemplate MQMulti; -PLUMED_REGISTER_ACTION(MQMulti,"MASS_VECTOR") -PLUMED_REGISTER_ACTION(MQMulti,"CHARGE_VECTOR") +typedef SelectMassCharge SelectMass; +typedef SelectMassCharge SelectCharge; +PLUMED_REGISTER_ACTION(SelectMass,"MASS_SCALAR") +PLUMED_REGISTER_ACTION(SelectCharge,"CHARGE_SCALAR") + +typedef ColvarShortcut MassShortcut; +typedef ColvarShortcut ChargeShortcut; +PLUMED_REGISTER_ACTION(MassShortcut,"MASS") +PLUMED_REGISTER_ACTION(ChargeShortcut,"CHARGE") + +typedef MultiColvarTemplate MassMulti; +typedef MultiColvarTemplate ChargeMulti; +PLUMED_REGISTER_ACTION(MassMulti,"MASS_VECTOR") +PLUMED_REGISTER_ACTION(ChargeMulti,"CHARGE_VECTOR") -void SelectMassCharge::registerKeywords( Keywords& keys ) { +template +void SelectMassCharge::registerKeywords( Keywords& keys ) { Colvar::registerKeywords( keys ); keys.add("atoms","ATOM","the atom number"); keys.add("atoms","ATOMS","the atom numbers that you would like to store the masses and charges of"); @@ -110,52 +125,89 @@ void SelectMassCharge::registerKeywords( Keywords& keys ) { keys.setDisplayName( acname.substr(0,und) ); keys.setValueDescription("the " + keys.getDisplayName() + " of the atom"); } -SelectMassCharge::SelectMassCharge(const ActionOptions&ao): +template +SelectMassCharge::SelectMassCharge(const ActionOptions&ao): PLUMED_COLVAR_INIT(ao) { - std::vector atoms; parseAtomList(-1,atoms,this); + std::vector atoms; + parseAtomList(-1,atoms,this); + theAtom=atoms[0]; /*Modetype mode=*/getModeAndSetupValues(this); requestAtoms(atoms); } -void SelectMassCharge::parseAtomList( int const num, std::vector& t, ActionAtomistic* aa ) { +template +void SelectMassCharge::parseAtomList( int const num, std::vector& t, ActionAtomistic* aa ) { aa->parseAtomList("ATOM",num,t); - if( t.size()==1 ) aa->log.printf(" for atom %d\n",t[0].serial()); - else if( num<0 || t.size()!=0 ) aa->error("Number of specified atoms should be 1"); + if( t.size()==1 ) { + aa->log.printf(" for atom %d\n",t[0].serial()); + } else if( num<0 || t.size()!=0 ) { + aa->error("Number of specified atoms should be 1"); + } } -SelectMassCharge::Modetype SelectMassCharge::getModeAndSetupValues( ActionWithValue* av ) { - av->addValueWithDerivatives(); av->setNotPeriodic(); bool constant=true; - ActionAtomistic* aa=dynamic_cast( av ); plumed_assert( aa ); +template +typename SelectMassCharge::Modetype SelectMassCharge::getModeAndSetupValues( ActionWithValue* av ) { + av->addValueWithDerivatives(); + av->setNotPeriodic(); + bool constant=true; + ActionAtomistic* aa=dynamic_cast( av ); + plumed_assert( aa ); for(unsigned i=0; igetNumberOfAtoms(); ++i) { std::pair p = aa->getValueIndices( aa->getAbsoluteIndex(i) ); - if( av->getName().find("MASS")!=std::string::npos && !aa->masv[p.first]->isConstant() ) constant=false; - if( av->getName().find("CHARGE")!=std::string::npos && !aa->chargev[p.first]->isConstant() ) constant=false; + if constexpr( mq == MC::Mass ) { + constant = aa->isMassConstant(p.first); + } else { + constant = aa->isChargeConstant(p.first); + } + } + if( !constant ) { + av->error("cannot deal with non-constant " + av->getName() + " values"); } - if( !constant ) av->error("cannot deal with non-constant " + av->getName() + " values"); (av->copyOutput(0))->setConstant(); return {}; } // calculator -void SelectMassCharge::calculate() { - std::vector masses(1), charges(1), value(1); - std::vector pos; std::vector > derivs; std::vector virial; - calculateCV( {}, masses, charges, pos, multiColvars::Ouput(value, derivs, virial), this ); setValue( value[0] ); + +template +void SelectMassCharge::calculate() { + std::vector posdummy; + std::vector > derivsdummy; + std::vector virialdummy; + + std::vector massesOrCharges(1); + if constexpr( mq == MC::Mass ) { + massesOrCharges[0]=getMass(theAtom.index()); + } else { + massesOrCharges[0]=getCharge(theAtom.index()); + } + std::vector vals(1); + + calculateCV( {}, massesOrCharges, massesOrCharges, posdummy, multiColvars::Ouput(vals, derivsdummy, virialdummy), this ); + + setValue( vals[0] ); + // does the code above give the same result as doing: + // if constexpr( mq == MC::Mass ) { + // setValue( getMass(theAtom.index) ); + // } else { + // setValue( getCharge(theAtom.index) ); + // } + // ??? + //calculateCV copies the first element from masses or charges into vals[0] } -void SelectMassCharge::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, - const std::vector& pos,multiColvars::Ouput out, const ActionAtomistic* aa ) { +template +void SelectMassCharge::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, + const std::vector& pos,multiColvars::Ouput out, const ActionAtomistic* aa ) { auto & vals=out.vals(); - if( aa->getName().find("MASSES")!=std::string::npos ) { + if constexpr(mq == MC::Mass) { vals[0]=masses[0]; - } else if( aa->chargesWereSet ) { + // } else if( aa->chargesWereSet ) { this is done in the getModeAndSetupValues by isChargeConstant + } else { vals[0]=charges[0]; } } -} -} - - - +} // namespace colvar +} // namespace PLMD diff --git a/src/core/ActionAtomistic.h b/src/core/ActionAtomistic.h index da8c1266b7..f04f08452e 100644 --- a/src/core/ActionAtomistic.h +++ b/src/core/ActionAtomistic.h @@ -35,10 +35,6 @@ namespace PLMD { class Pbc; class PDB; -namespace colvar { -class SelectMassCharge; -} - /// \ingroup MULTIINHERIT /// Action used to create objects that access the positions of the atoms from the MD code class ActionAtomistic : @@ -46,7 +42,6 @@ class ActionAtomistic : { friend class Group; friend class DomainDecomposition; - friend class colvar::SelectMassCharge; friend class ActionWithVirtualAtom; std::vector indexes; // the set of needed atoms @@ -82,9 +77,9 @@ class ActionAtomistic : protected: bool chargesWereSet; void setExtraCV(const std::string &name); +public: /// Used to interpret whether this index is a virtual atom or a real atom std::pair getValueIndices( const AtomNumber& i ) const ; -public: /// Request an array of atoms. /// This method is used to ask for a list of atoms. Atoms /// should be asked for by number. If this routine is called @@ -116,8 +111,12 @@ class ActionAtomistic : const double & getEnergy()const; /// Get mass of i-th atom double getMass(int i)const; +/// Check if the mass of the i-th atom is constant + bool isMassConstant(int i) const; /// Get charge of i-th atom double getCharge(int i)const; +/// Check if the charge of the i-th atom is constant + bool isChargeConstant(int i) const; /// Get the force acting on a particular atom Vector getForce( const std::pair& a ) const ; /// Add force to an atom @@ -205,12 +204,23 @@ double ActionAtomistic::getMass(int i)const { return masses[i]; } +inline +bool ActionAtomistic::isMassConstant(int i)const { + return masv[i]->isConstant(); +} + inline double ActionAtomistic::getCharge(int i) const { if( !chargesWereSet ) error("charges were not passed to plumed"); return charges[i]; } +inline +bool ActionAtomistic::isChargeConstant(int i) const { + if( !chargesWereSet ) error("charges were not passed to plumed"); + return chargev[i]->isConstant(); +} + inline const std::vector & ActionAtomistic::getAbsoluteIndexes()const { return indexes; From 3918864a08b00b5c8ed2c7f619ba5e527bd14f1d Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 20 Nov 2024 13:10:03 +0100 Subject: [PATCH 07/10] forgot to remove a duplicated declaration in quaternion --- src/crystdistrib/Quaternion.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/crystdistrib/Quaternion.cpp b/src/crystdistrib/Quaternion.cpp index 86fe8ac21a..d6f05bc327 100644 --- a/src/crystdistrib/Quaternion.cpp +++ b/src/crystdistrib/Quaternion.cpp @@ -103,8 +103,6 @@ class Quaternion : public Colvar { public: static void registerKeywords( Keywords& keys ); explicit Quaternion(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - // active methods: void calculate() override; MULTICOLVAR_DEFAULT(::PLMD::colvar::multiColvars::emptyMode); From 060b5d3541c7840fba6b8175db7845193be5910a Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 20 Nov 2024 14:37:03 +0100 Subject: [PATCH 08/10] removing some ordering mishap --- src/colvar/Torsion.cpp | 4 ++-- src/crystdistrib/Quaternion.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp index 2a9c6fb314..36b249f91a 100644 --- a/src/colvar/Torsion.cpp +++ b/src/colvar/Torsion.cpp @@ -145,10 +145,10 @@ void Torsion::registerKeywords(Keywords& keys) { Torsion::Torsion(const ActionOptions&ao): PLUMED_COLVAR_INIT(ao), - pbc(true), value(1), derivs(1), - virial(1) + virial(1), + pbc(true) { derivs[0].resize(6); std::vector atoms; std::vector v1; ActionAtomistic::parseAtomList("VECTOR1",v1); diff --git a/src/crystdistrib/Quaternion.cpp b/src/crystdistrib/Quaternion.cpp index d6f05bc327..8a3edc600a 100644 --- a/src/crystdistrib/Quaternion.cpp +++ b/src/crystdistrib/Quaternion.cpp @@ -145,7 +145,7 @@ Quaternion::Quaternion(const ActionOptions&ao): requestAtoms(atoms); } -void Quaternion::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { +void Quaternion::parseAtomList(int num, std::vector& t, ActionAtomistic* aa ) { aa->parseAtomList("ATOMS",num,t); if( t.size()==3 ) aa->log.printf(" involving atoms %d %d %d\n",t[0].serial(),t[1].serial(),t[0].serial()); } From 0ef2ad0ef5c41c30790c7f6658f58163c323f708 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti Date: Thu, 21 Nov 2024 08:47:21 +0100 Subject: [PATCH 09/10] setting up Input setting up a (very complex) input variable adding a builder pattern feature, that may be useful runned atyle adding tests for Input set tests for baisc input nearly everithing tested simplifying the Input interface --- regtest/multicolvar/rt-input/Makefile | 1 + regtest/multicolvar/rt-input/config | 1 + regtest/multicolvar/rt-input/main.cpp | 201 ++++++++++++++++++ regtest/multicolvar/rt-input/output.reference | 66 ++++++ src/colvar/Angle.cpp | 12 +- src/colvar/DihedralCorrelation.cpp | 9 +- src/colvar/Dipole.cpp | 24 +-- src/colvar/Distance.cpp | 9 +- src/colvar/MultiColvarTemplate.h | 134 ++++++++++-- src/colvar/Plane.cpp | 10 +- src/colvar/Position.cpp | 12 +- src/colvar/SelectMassCharge.cpp | 22 +- src/colvar/Torsion.cpp | 10 +- src/crystdistrib/Quaternion.cpp | 9 +- 14 files changed, 449 insertions(+), 71 deletions(-) create mode 100644 regtest/multicolvar/rt-input/Makefile create mode 100644 regtest/multicolvar/rt-input/config create mode 100644 regtest/multicolvar/rt-input/main.cpp create mode 100644 regtest/multicolvar/rt-input/output.reference diff --git a/regtest/multicolvar/rt-input/Makefile b/regtest/multicolvar/rt-input/Makefile new file mode 100644 index 0000000000..3703b27cea --- /dev/null +++ b/regtest/multicolvar/rt-input/Makefile @@ -0,0 +1 @@ +include ../../scripts/test.make diff --git a/regtest/multicolvar/rt-input/config b/regtest/multicolvar/rt-input/config new file mode 100644 index 0000000000..df1f95bf3e --- /dev/null +++ b/regtest/multicolvar/rt-input/config @@ -0,0 +1 @@ +type=make diff --git a/regtest/multicolvar/rt-input/main.cpp b/regtest/multicolvar/rt-input/main.cpp new file mode 100644 index 0000000000..d0dd3fdc5e --- /dev/null +++ b/regtest/multicolvar/rt-input/main.cpp @@ -0,0 +1,201 @@ +#include "plumed/colvar/MultiColvarTemplate.h" +#include "plumed/tools/Vector.h" + +#include +#include +#include +#include + +using namespace PLMD; +using colvar::multiColvars::Input; + +template void pv(std::ofstream &stream, const std::vector &v) { + // todo: put this is a common header + std::string div = ""; + stream << "(" << v.size() << ") "; + for (const auto &t : v) { + stream << div << t; + div = ", "; + } +} + +void print(std::ofstream &stream, const std::vector &masses, + const std::vector &charges, + const std::vector positions) { + stream << "Masses: "; + pv(stream, masses); + stream << "\nCharges: "; + pv(stream, charges); + stream << "\nPositions: "; + pv(stream, positions); + stream << "\n"; +} + +void inputTest(std::ofstream &stream, Input input) { + try { + const auto &masses = input.masses(); + stream << " Masses: "; + pv(stream, masses); + stream << "\n"; + } catch (std::bad_variant_access &) { + stream << " There are no masses\n"; + } + + try { + const auto &charges = input.charges(); + stream << " Charges: "; + pv(stream, charges); + stream << "\n"; + } catch (std::bad_variant_access &) { + stream << " There are no charges\n"; + } + try { + const auto &positions = input.positions(); + stream << " Positions: "; + pv(stream, positions); + stream << "\n"; + } catch (std::bad_variant_access &) { + stream << " There are no positions\n"; + } + stream << std::endl; +} + +void workOnMass(Input input) { + auto &masses = input.var_masses(); + for (auto &m : masses) { + m *= 2; + } +} + +void workOnCharges(Input input) { + auto &charges = input.var_charges(); + for (auto &c : charges) { + c *= 2; + } +} + +void workOnPositions(Input input) { + auto &pos = input.var_positions(); + for (auto &p : pos) { + p = p * 2.0; + } +} + +class testData { + using vd = std::vector; + using vv = std::vector; + vv p = {{0, 0, 0}, {1, 1, 1}}; + vd c = {-1, 1}; + vd m = {3, 4}; + +public: + testData() = default; + vd &masses() { return m; } + vd &charges() { return c; } + vv &positions() { return p; } + const vd &c_masses() const { return m; } + const vd &c_charges() const { return c; } + const vv &c_positions() const { return p; } +}; + +class test { + using vd = std::vector; + using vv = std::vector; + +public: + test() = default; + void dotestWithReference(std::ofstream &output) { + testData td; + output << "*dotestWithReference:" << std::endl; + output << "*Passing full input:\n"; + inputTest(output, Input(td.masses(), td.charges(), td.positions())); + output << "*Passing only positions:\n"; + inputTest(output, Input().positions(td.positions())); + output << "*Passing only masses:\n"; + inputTest(output, Input().masses(td.masses())); + output << "*Passing only charges:\n"; + inputTest(output, Input().charges(td.charges())); + } + void dotestWithConstReference(std::ofstream &output) { + testData td; + output << "*dotestWithConstReference:" << std::endl; + output << "*Passing full input:\n"; + inputTest(output, Input(td.c_masses(), td.c_charges(), td.c_positions())); + output << "*Passing only positions:\n"; + inputTest(output, Input().positions(td.c_positions())); + output << "*Passing only masses:\n"; + inputTest(output, Input().masses(td.c_masses())); + output << "*Passing only charges:\n"; + inputTest(output, Input().charges(td.c_charges())); + } + + void testUsingInputAsBuffer_masses(std::ofstream &output) { + testData td; + output << "*testUsingInputAsBuffer_masses:" << std::endl; + output << "*Passing constant_reference:\n"; + try { + workOnMass(Input(td.c_masses(), td.c_charges(), td.c_positions())); + } catch (std::bad_variant_access &) { + output << " should throw\n"; + } + output << "*Passing reference:\n"; + output << " Initial values\n Masses:"; + pv(output, td.masses()); + workOnMass(Input(td.masses(), td.charges(), td.positions())); + output << "\n"; + // this extra line is to remember that something has changed!!! + output << " Values after call\n Masses:"; + pv(output, td.masses()); + output << std::endl; + } + void testUsingInputAsBuffer_charges(std::ofstream &output) { + testData td; + output << "*testUsingInputAsBuffer_charges:" << std::endl; + output << "*Passing constant_reference:\n"; + try { + workOnCharges(Input(td.c_masses(), td.c_charges(), td.c_positions())); + } catch (std::bad_variant_access &) { + output << " should throw\n"; + } + output << "*Passing reference:\n"; + output << " Initial values\n Charges:"; + pv(output, td.c_charges()); + workOnCharges(Input(td.masses(), td.charges(), td.positions())); + output << "\n"; + // this extra line is to remember that something has changed!!! + output << " Values after call\n Charges:"; + pv(output, td.c_charges()); + output << std::endl; + } + void testUsingInputAsBuffer_positions(std::ofstream &output) { + testData td; + output << "*testUsingInputAsBuffer_positions:" << std::endl; + output << "*Passing constant_reference:\n"; + try { + workOnPositions(Input(td.c_masses(), td.c_charges(), td.c_positions())); + } catch (std::bad_variant_access &) { + output << " should throw\n"; + } + output << "*Passing reference:\n"; + output << " Initial values\n Positions:"; + pv(output, td.positions()); + workOnPositions(Input(td.masses(), td.charges(), td.positions())); + output << "\n"; + // this extra line is to remember that something has changed!!! + output << " Values after call\n Positions:"; + pv(output, td.positions()); + output << std::endl; + } +}; + +int main() { + std::ofstream output("output"); + + test t; + t.dotestWithReference(output); + t.dotestWithConstReference(output); + t.testUsingInputAsBuffer_masses(output); + t.testUsingInputAsBuffer_charges(output); + t.testUsingInputAsBuffer_positions(output); + return 0; +} diff --git a/regtest/multicolvar/rt-input/output.reference b/regtest/multicolvar/rt-input/output.reference new file mode 100644 index 0000000000..43c6532d27 --- /dev/null +++ b/regtest/multicolvar/rt-input/output.reference @@ -0,0 +1,66 @@ +*dotestWithReference: +*Passing full input: + Masses: (2) 3, 4 + Charges: (2) -1, 1 + Positions: (2) 0 0 0, 1 1 1 + +*Passing only positions: + There are no masses + There are no charges + Positions: (2) 0 0 0, 1 1 1 + +*Passing only masses: + Masses: (2) 3, 4 + There are no charges + There are no positions + +*Passing only charges: + There are no masses + Charges: (2) -1, 1 + There are no positions + +*dotestWithConstReference: +*Passing full input: + Masses: (2) 3, 4 + Charges: (2) -1, 1 + Positions: (2) 0 0 0, 1 1 1 + +*Passing only positions: + There are no masses + There are no charges + Positions: (2) 0 0 0, 1 1 1 + +*Passing only masses: + Masses: (2) 3, 4 + There are no charges + There are no positions + +*Passing only charges: + There are no masses + Charges: (2) -1, 1 + There are no positions + +*testUsingInputAsBuffer_masses: +*Passing constant_reference: + should throw +*Passing reference: + Initial values + Masses:(2) 3, 4 + Values after call + Masses:(2) 6, 8 +*testUsingInputAsBuffer_charges: +*Passing constant_reference: + should throw +*Passing reference: + Initial values + Charges:(2) -1, 1 + Values after call + Charges:(2) -2, 2 +*testUsingInputAsBuffer_positions: +*Passing constant_reference: + should throw +*Passing reference: + Initial values + Positions:(2) 0 0 0, 1 1 1 + Values after call + Positions:(2) 0 0 0, 2 2 2 diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp index f5b9e4dd38..7f458330df 100644 --- a/src/colvar/Angle.cpp +++ b/src/colvar/Angle.cpp @@ -100,7 +100,7 @@ Calculate multiple angles. class Angle : public Colvar { bool pbc; - std::vector value, masses, charges; + std::vector value; std::vector > derivs; std::vector virial; public: @@ -168,19 +168,21 @@ Angle::Angle(const ActionOptions&ao): void Angle::calculate() { if(pbc) makeWhole(); - calculateCV( {}, masses, charges, getPositions(), multiColvars::Ouput(value, derivs, virial), this ); + calculateCV( {}, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this ); setValue( value[0] ); for(unsigned i=0; i& /*masses*/, const std::vector& /*charges*/, - const std::vector& pos, - multiColvars::Ouput out, const ActionAtomistic* aa ) { +void Angle::calculateCV( Modetype /*mode*/, + multiColvars::Input const in, + multiColvars::Ouput out, + const ActionAtomistic* aa ) { auto & vals=out.vals(); auto & derivs=out.derivs(); auto & virial=out.virial(); + const auto & pos=in.positions(); Vector dij,dik; dij=delta(pos[2],pos[3]); dik=delta(pos[1],pos[0]); diff --git a/src/colvar/DihedralCorrelation.cpp b/src/colvar/DihedralCorrelation.cpp index 1a5e42c604..41ac7388bd 100644 --- a/src/colvar/DihedralCorrelation.cpp +++ b/src/colvar/DihedralCorrelation.cpp @@ -64,7 +64,7 @@ Measure the correlation between a multiple pairs of dihedral angles class DihedralCorrelation : public Colvar { private: bool pbc; - std::vector value, masses, charges; + std::vector value; std::vector > derivs; std::vector virial; public: @@ -124,18 +124,19 @@ DihedralCorrelation::Modetype DihedralCorrelation::getModeAndSetupValues( Action void DihedralCorrelation::calculate() { if(pbc) makeWhole(); - calculateCV( {}, masses, charges, getPositions(), multiColvars::Ouput(value, derivs, virial), this ); + calculateCV( {}, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this ); setValue( value[0] ); for(unsigned i=0; i& masses, const std::vector& charges, - const std::vector& pos, +void DihedralCorrelation::calculateCV(Modetype /*mode*/, + multiColvars::Input const in, multiColvars::Ouput out, const ActionAtomistic* aa ) { auto & vals=out.vals(); auto & derivs=out.derivs(); auto & virial=out.virial(); + const auto& pos = in.positions(); const Vector d10=delta(pos[1],pos[0]); const Vector d11=delta(pos[2],pos[1]); const Vector d12=delta(pos[3],pos[2]); diff --git a/src/colvar/Dipole.cpp b/src/colvar/Dipole.cpp index 972fc92b33..2ef27f7736 100644 --- a/src/colvar/Dipole.cpp +++ b/src/colvar/Dipole.cpp @@ -83,17 +83,11 @@ class Dipole : public Colvar { explicit Dipole(const ActionOptions&); void calculate() override; static void registerKeywords(Keywords& keys); - MULTICOLVAR_SETTINGS_BASE(multiColvars::components); - //Declaring this playinly becasue we are using modifying the charges - static void calculateCV( Modetype mode, - const std::vector& masses, - std::vector& charges, - const std::vector& pos, - multiColvars::Ouput out, - const ActionAtomistic* aa ); + MULTICOLVAR_DEFAULT(multiColvars::components); private: std::vector ga_lista; - std::vector value, masses, charges; + std::vector value; + std::vector charges; std::vector > derivs; std::vector virial; Value* valuex=nullptr; @@ -179,7 +173,9 @@ void Dipole::calculate() if(!nopbc) makeWhole(); unsigned N=getNumberOfAtoms(); for(unsigned i=0; i& masses, std::vector& charges, - const std::vector& pos,multiColvars::Ouput out, const ActionAtomistic* aa ) { +void Dipole::calculateCV( Modetype mode, + multiColvars::Input const in, + multiColvars::Ouput out, + const ActionAtomistic* aa ) { auto & vals=out.vals(); auto & derivs=out.derivs(); auto & virial=out.virial(); + const auto & pos=in.positions(); + auto & charges=in.var_charges(); unsigned N=pos.size(); double ctot=0.; for(unsigned i=0; i value{{0.0}}; - std::vectormasses{{0.0}}; - std::vector charges{{0.0}}; std::vector > derivs{std::vector{Vector{},Vector{}}}; std::vector virial{Tensor()}; Modetype mode_; @@ -221,7 +219,7 @@ Distance::Modetype Distance::getModeAndSetupValues( ActionWithValue* av ) { void Distance::calculate() { if(pbc) makeWhole(); - calculateCV( mode_, masses, charges, getPositions(), multiColvars::Ouput(value, derivs, virial), this ); + calculateCV( mode_, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this ); switch (mode_) { case Modetype::withCompontents: { @@ -278,12 +276,13 @@ void Distance::calculate() { } } -void Distance::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, +void Distance::calculateCV( Modetype mode, + multiColvars::Input const in, multiColvars::Ouput out, const ActionAtomistic* aa ) { auto & vals=out.vals(); auto & derivs=out.derivs(); auto & virial=out.virial(); + const auto & pos = in.positions(); Vector distance=delta(pos[0],pos[1]); const double value=distance.modulo(); const double invvalue=1.0/value; diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h index fb5d6c386c..fb1dbcbad6 100644 --- a/src/colvar/MultiColvarTemplate.h +++ b/src/colvar/MultiColvarTemplate.h @@ -22,6 +22,9 @@ #ifndef __PLUMED_colvar_MultiColvarTemplate_h #define __PLUMED_colvar_MultiColvarTemplate_h +#include +#include + #include "core/ActionWithVector.h" namespace PLMD { @@ -45,12 +48,10 @@ enum class scaledComponents { //TODO: test this -class Ouput -{ - //these are const pointers to non const data,meaning YOU can't change the pointer - std::vector * vals_=nullptr; +class Ouput final { + std::vector* vals_=nullptr; std::vector >* derivs_=nullptr; - std::vector * virial_=nullptr; + std::vector* virial_=nullptr; Ouput()=default; public: //const because the data in the class is not changing: we are modifying data pointed by the const pointers @@ -60,9 +61,15 @@ class Ouput Ouput(std::vector& vals, std::vector >& derivs, std::vector& virial) - : vals_(&vals), derivs_(&derivs), virial_(&virial) {} - Ouput(Ouput const& other) : vals_(other.vals_), derivs_(other.derivs_), virial_(other.virial_) {}; - Ouput(Ouput && other) noexcept : Ouput() { + : vals_(&vals), + derivs_(&derivs), + virial_(&virial) {} + Ouput(Ouput const& other) + : vals_(other.vals_), + derivs_(other.derivs_), + virial_(other.virial_) {}; + Ouput(Ouput && other) noexcept + : Ouput() { swap(*this, other); }; // the input is initialized as copy or move ctor @@ -80,6 +87,104 @@ class Ouput ~Ouput() {}//this does not own anything }; +class Input final { + using vd=std::vector*; + using cvd = const std::vector*; + using vv=std::vector*; + using cvv = const std::vector*; + std::variant*,const std::vector*> masses_; + std::variant*,const std::vector*> charges_; + std::variant*,const std::vector*> positions_; +public: +//const because the data in the class is not changing: we are modifying data pointed by the const pointers +//get_if return nullptr on error +//get throws on error +//references to variable data, work only if the data passed is a reference to variable data + std::vector&var_masses() const {return *std::get*>(masses_);} + std::vector&var_charges() const {return *std::get*>(charges_);} + std::vector&var_positions() const {return *std::get*>(positions_);} + //references to constant data, must ALWAYS work + const std::vector&masses() const { + if(auto *p=std::get_if(&masses_)) { + return **p; + } + //get_if takes a ptr, get takes a ref... + //this throws if there is no data + return *std::get(masses_); + } + const std::vector&charges() const { + if(auto *p=std::get_if(&charges_)) { + return **p; + } + //get_if takes a ptr, get takes a ref... + //this throws if there is no data + return *std::get(charges_); + } + const std::vector&positions() const { + if(auto *p=std::get_if(&positions_)) { + return **p; + } + //get_if takes a ptr, get takes a ref... + //this throws if there is no data + return *std::get(positions_); + } + + Input()=default; + template + Input(masses_T& masses, + charges_T& charges, + positions_T& positions) + : masses_(&masses), + charges_(&charges), + positions_(&positions) { + static_assert(std::is_same_v> || std::is_same_v>); + static_assert(std::is_same_v> || std::is_same_v>); + static_assert(std::is_same_v> || std::is_same_v>); + } + Input(Input const& other) + : masses_(other.masses_), + charges_(other.charges_), + positions_(other.positions_) {}; + Input(Input && other) noexcept + : Input() { + swap(*this, other); + }; + //with the builder pattern I ma not need the ubercomplex ctor and the just things that I set up... + Input& charges(std::vector& charges) { + charges_=&charges; + return *this; + } + Input& charges(const std::vector& charges) { + charges_=&charges; + return *this; + } + Input& masses(std::vector& masses) { + masses_=&masses; + return *this; + } + Input& masses(const std::vector& masses) { + masses_=&masses; + return *this; + } + Input& positions(std::vector& positions) { + positions_=&positions; + return *this; + } + Input& positions(const std::vector& positions) { + positions_=&positions; + return *this; + } + + // the input is initialized as copy or move ctor + Input& operator=(Input other) =delete; + friend void swap(Input& a, Input& b) noexcept { + std::swap(a.masses_, b.masses_); + std::swap(a.charges_, b.charges_); + std::swap(a.positions_, b.positions_); + } + ~Input() {}//this does not own anything +}; + } // namespace multiColvars /*MANUAL DRAFT @@ -97,7 +202,9 @@ To setup a CV compatible with multicolvartemplate you need to add - these functions will need to match the foloowing signatures (withour the constness of the non const& arguments): - `static void parseAtomList( int num, std::vector& t, ActionAtomistic* aa );` - `Modetype getModeAndSetupValues( ActionWithValue* av )` - - `static void calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, const std::vector& pos, PLMD::colvars::multiColvars::Ouput out, const ActionAtomistic* aa ) + - `static void calculateCV( Modetype mode, PLMD::colvar::multiColvars::Input in, PLMD::colvar::multiColvars::Ouput out, const ActionAtomistic* aa ) + - the input variable contain the references to position, masses, charges + - the output variable contain the references to value, derivs, virial,and acts as a return argument - this function will be called by MultiColvarTemplate to calculate the CV value on the inputs - to avoid code repetitition you should change ::calculate() to call this function - By default all the inputs are const ref, but the constedness can be changed, since the MulticolvarTemplate will pass the plain references @@ -110,17 +217,12 @@ To setup a CV compatible with multicolvartemplate you need to add static void parseAtomList( int num, std::vector& t, ActionAtomistic* aa ); \ static Modetype getModeAndSetupValues( ActionWithValue* av ) #define MULTICOLVAR_SETTINGS_CALCULATE_CONST() static void calculateCV( Modetype mode, \ - const std::vector& masses, \ - const std::vector& charges, \ - const std::vector& pos, \ + PLMD::colvar::multiColvars::Input in, \ PLMD::colvar::multiColvars::Ouput out, \ const ActionAtomistic* aa ) #define MULTICOLVAR_DEFAULT(type) MULTICOLVAR_SETTINGS_BASE(type); \ MULTICOLVAR_SETTINGS_CALCULATE_CONST() - - - //^no ';' here so that the macro will "consume" it and not generate a double semicolon error template @@ -285,7 +387,7 @@ void MultiColvarTemplate::performTask( const unsigned& task_index, MultiValu } } // Calculate the CVs using the method in the Colvar - CV::calculateCV( mode, mass, charge, fpositions, multiColvars::Ouput{values, derivs, virial}, this ); + CV::calculateCV( mode, multiColvars::Input(mass, charge, fpositions), multiColvars::Ouput{values, derivs, virial}, this ); for(unsigned i=0; i value, masses, charges; + std::vector value; std::vector > derivs; std::vector virial; public: @@ -131,17 +131,19 @@ Plane::Plane(const ActionOptions&ao): void Plane::calculate() { if(pbc) makeWhole(); - calculateCV( {}, masses, charges, getPositions(), multiColvars::Ouput(value, derivs, virial), this ); + calculateCV( {}, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this ); setValue( value[0] ); for(unsigned i=0; i& masses, const std::vector& charges, - const std::vector& pos,multiColvars::Ouput out, const ActionAtomistic* aa ) { +void Plane::calculateCV( Modetype /*mode*/, + multiColvars::Input const in, + multiColvars::Ouput out, const ActionAtomistic* aa ) { auto & vals=out.vals(); auto & derivs=out.derivs(); auto & virial=out.virial(); + const auto & pos = in.positions(); Vector d1=delta( pos[1], pos[0] ); Vector d2=delta( pos[2], pos[3] ); Vector cp = crossProduct( d1, d2 ); diff --git a/src/colvar/Position.cpp b/src/colvar/Position.cpp index 9e902befd1..0abf9c016e 100644 --- a/src/colvar/Position.cpp +++ b/src/colvar/Position.cpp @@ -101,8 +101,6 @@ class Position : public Colvar { Modetype components; bool pbc; std::vector value; - std::vector masses; - std::vector charges; std::vector > derivs; std::vector virial; }; @@ -181,7 +179,7 @@ void Position::calculate() { } else { distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0)); } - calculateCV( components, masses, charges, distance, multiColvars::Ouput(value, derivs, virial), this ); + calculateCV( components, multiColvars::Input().positions(distance), multiColvars::Ouput(value, derivs, virial), this ); switch (components) { case Modetype::scaled: { Value* valuea=getPntrToComponent("a"); @@ -216,12 +214,16 @@ void Position::calculate() { } } -void Position::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, +void Position::calculateCV( Modetype mode, + multiColvars::Input const in, multiColvars::Ouput out, const ActionAtomistic* aa ) { auto & vals=out.vals(); auto & derivs=out.derivs(); auto & virial=out.virial(); + const auto &pos = in.positions(); + //////////////////////////////////////////////////////////////////////////// + //why here there is not the distance from 000 treatment?????? + /////////////////////////////////////////////////////////////////////////// switch (mode) { case Modetype::scaled: { diff --git a/src/colvar/SelectMassCharge.cpp b/src/colvar/SelectMassCharge.cpp index e396a4c646..69a80d4fd9 100644 --- a/src/colvar/SelectMassCharge.cpp +++ b/src/colvar/SelectMassCharge.cpp @@ -172,19 +172,18 @@ typename SelectMassCharge::Modetype SelectMassCharge::getModeAndSetupVal template void SelectMassCharge::calculate() { - std::vector posdummy; + std::vector > derivsdummy; std::vector virialdummy; + std::vector vals(1); - std::vector massesOrCharges(1); if constexpr( mq == MC::Mass ) { - massesOrCharges[0]=getMass(theAtom.index()); + std::vector mass{getMass(theAtom.index())}; + calculateCV( {}, multiColvars::Input().masses(mass), multiColvars::Ouput(vals, derivsdummy, virialdummy), this ); } else { - massesOrCharges[0]=getCharge(theAtom.index()); + std::vector charge{getCharge(theAtom.index())}; + calculateCV( {}, multiColvars::Input().charges(charge), multiColvars::Ouput(vals, derivsdummy, virialdummy), this ); } - std::vector vals(1); - - calculateCV( {}, massesOrCharges, massesOrCharges, posdummy, multiColvars::Ouput(vals, derivsdummy, virialdummy), this ); setValue( vals[0] ); // does the code above give the same result as doing: @@ -198,14 +197,15 @@ void SelectMassCharge::calculate() { } template -void SelectMassCharge::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, - const std::vector& pos,multiColvars::Ouput out, const ActionAtomistic* aa ) { +void SelectMassCharge::calculateCV( Modetype /*mode*/, + multiColvars::Input const in, + multiColvars::Ouput out, const ActionAtomistic* aa ) { auto & vals=out.vals(); if constexpr(mq == MC::Mass) { - vals[0]=masses[0]; + vals[0]= in.masses()[0]; // } else if( aa->chargesWereSet ) { this is done in the getModeAndSetupValues by isChargeConstant } else { - vals[0]=charges[0]; + vals[0]=in.charges()[0]; } } diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp index 36b249f91a..83205ba674 100644 --- a/src/colvar/Torsion.cpp +++ b/src/colvar/Torsion.cpp @@ -117,7 +117,7 @@ class Torsion : public Colvar { }; MULTICOLVAR_DEFAULT(torsionModes); private: - std::vector value, masses, charges; + std::vector value; std::vector > derivs; std::vector virial; bool pbc; @@ -228,7 +228,7 @@ void Torsion::calculate() { if(pbc) { makeWhole(); } - calculateCV( mode, masses, charges, getPositions(), multiColvars::Ouput(value, derivs, virial), this ); + calculateCV( mode, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this ); for(unsigned i=0; i<6; ++i) { setAtomsDerivatives(i,derivs[0][i] ); } @@ -236,13 +236,13 @@ void Torsion::calculate() { setBoxDerivatives( virial[0] ); } -void Torsion::calculateCV( Modetype mode, const std::vector& masses, const std::vector& charges, - const std::vector& pos, +void Torsion::calculateCV( Modetype mode, + multiColvars::Input const in, multiColvars::Ouput out, const ActionAtomistic* aa ) { auto & vals=out.vals(); auto & derivs=out.derivs(); auto & virial=out.virial(); - + const auto & pos = in.positions(); Vector d0=delta(pos[1],pos[0]); Vector d1=delta(pos[3],pos[2]); Vector d2=delta(pos[5],pos[4]); diff --git a/src/crystdistrib/Quaternion.cpp b/src/crystdistrib/Quaternion.cpp index 8a3edc600a..02975cecba 100644 --- a/src/crystdistrib/Quaternion.cpp +++ b/src/crystdistrib/Quaternion.cpp @@ -97,7 +97,7 @@ See \ref QUATERNION for more details class Quaternion : public Colvar { private: bool pbc; - std::vector value, masses, charges; + std::vector value; std::vector > derivs; std::vector virial; public: @@ -162,7 +162,7 @@ Quaternion::Modetype Quaternion::getModeAndSetupValues( ActionWithValue* av ) { void Quaternion::calculate() { if(pbc) makeWhole(); - calculateCV( {}, masses, charges, getPositions(), PLMD::colvar::multiColvars::Ouput(value, derivs, virial), this ); + calculateCV( {}, PLMD::colvar::multiColvars::Input().positions(getPositions()), PLMD::colvar::multiColvars::Ouput(value, derivs, virial), this ); for(unsigned j=0; j<4; ++j) { Value* valuej=getPntrToComponent(j); for(unsigned i=0; i<3; ++i) setAtomsDerivatives(valuej,i,derivs[j][i] ); @@ -172,12 +172,13 @@ void Quaternion::calculate() { } // calculator -void Quaternion::calculateCV( Modetype /*mode*/, const std::vector& masses, const std::vector& charges, - const std::vector& pos, +void Quaternion::calculateCV( Modetype /*mode*/, + PLMD::colvar::multiColvars::Input const in, PLMD::colvar::multiColvars::Ouput out, const ActionAtomistic* aa ) { auto & vals=out.vals(); auto & derivs=out.derivs(); auto & virial=out.virial(); + const auto & pos = in.positions(); //declarations Vector vec1_comp = delta( pos[0], pos[1] ); //components between atom 1 and 2 Vector vec2_comp = delta( pos[0], pos[2] ); //components between atom 1 and 3 From 0f524d336cce6e6075ad450165733705eff1d813 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti Date: Fri, 22 Nov 2024 09:55:52 +0100 Subject: [PATCH 10/10] udating some dependencies --- regtest/crystdistrib/rt-dops-shortcut/config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/regtest/crystdistrib/rt-dops-shortcut/config b/regtest/crystdistrib/rt-dops-shortcut/config index 557dccbe32..d09ec688c6 100644 --- a/regtest/crystdistrib/rt-dops-shortcut/config +++ b/regtest/crystdistrib/rt-dops-shortcut/config @@ -1,4 +1,4 @@ type=driver -plumed_modules=crystdistrib +plumed_modules="crystdistrib adjmat" arg="--plumed plumed.dat --mf_xtc em.xtc --mc mcfile --timestep=0.001 --trajectory-stride=25" extra_files="../../trajectories/em.xtc"