plumed2 icon indicating copy to clipboard operation
plumed2 copied to clipboard

Help wanted: How to reference the Q6 and it's derivatives defined in Plumed?

Open tuoping opened this issue 3 years ago • 6 comments

Hello. I want to use the Q6 defined in Plumed. But failed to find it in wrapper. I have two questions. 1. Is it compiled in libplumed.so or libplumedkernel.so? If it is, how do I reference it? If it is not, is it possible for me to include it in the compiled static lib? Or is there a way for me to use it directly? Thank you sincerely.

tuoping avatar Jul 06 '21 07:07 tuoping

@tuoping it is not clear what you want to do. Do you want to call these functions from a C/C++ program?

  • The Q6 implementation is included in the crystallization module, that is not enabled by default (configure with --enable-modules=crystallization).
  • Once it's enabled, the resulting binary code will end up in libplumedKernel.so.
  • I am not sure there is a way to use the relevant functions from a C/C++ program. However, you should be able to compute the Q6 using plumed driver on an existing trajectory of to call PLUMED through its cmd interface (from C/C++/Python/Fortran) to compute any collective variable, including this one.

GiovanniBussi avatar Jul 06 '21 12:07 GiovanniBussi

Thank you for your kind reply. I do want to call the functions from a C++ program. And I found that every function requires a ActionOption parameter as input, which makes it hard to use. Do you think it's possible or usefull to detach Action related parts in these functions, so that these functions can be called independently?

tuoping avatar Jul 07 '21 01:07 tuoping

You are right. In some cases we have tools that can be used easily from outside PLUMED (e.g. RMSD), but most CVs are deeply integrated in our framework and difficult to use as standalone.

I would suggest that you create a PLMD::Plumed instance and use the PLMD::Plumed::cmd interface (the same one that is used from MD codes such as LAMMPS and NAMD) to communicate with PLUMED the coordinates on which you want to compute the Q6 and retrieve the CVs. There is an example in Python here, but the interface in Python and C++ is virtually identical. You can find many examples of how to use the C++ interface in our regtests with cd regtests; git grep "cmd(".

I recommend to use the master branch since there the C++ interface implements type checks (#653) and so it is much less likely to make mistakes. There is one specific tricky point, namely the getDataRank command was meant to be used from Python and you should thus feed it with long int* (and not int*).

Given the fact that variables such as Q6 are very expensive, I don't expect a significant overhead with respect to a hypothetical direct use of the Q6 class.

GiovanniBussi avatar Jul 07 '21 07:07 GiovanniBussi

This is very helpfull. Thank you so much.

tuoping avatar Jul 08 '21 05:07 tuoping

@tuoping please let us know if you succeeded (we are always interested in feedbacks from people using PLUMED differently from what we expected)

GiovanniBussi avatar Jul 11 '21 15:07 GiovanniBussi

Hello, @GiovanniBussi. Thank you for your help earlier this year. I managed to use MultiColvar, Q6 to be specific, by p->cmd(). And I did it by #include MultiColvarBase.h in PlumedMain.cpp, and build StoreDataStash in just_calculate().

void PlumedMain::justCalculate() {
...
multicolvar::MultiColvarBase* amc=dynamic_cast<multicolvar::MultiColvarBase*>(p);
if(amc){
    // log<<"> getNumberOfVessels(): "<<amc->getNumberOfVessels()<<"\n";
    // log<<"> nactive_tasks: "<<amc->getCurrentNumberOfActiveTasks()<<"\n";
    vesselbase::StoreDataVessel* mystash = amc->buildDataStashes( NULL );
    NumberOfTasks=amc->getFullNumberOfTasks();
    NumberOfQuantities=amc->getNumberOfQuantities();
    int NumberOfStoredValues=mystash->getNumberOfStoredValues();
    // std::vector<double>  myvalues( NumberOfQuantities );
    myvalues.resize(NumberOfTasks);
    log<<">> getNumberOfQuantities(): "<<amc->getNumberOfQuantities()<<"\n";
    log<<">> getNumberOfStoredValues(): "<<NumberOfStoredValues<<"\n";
    for(unsigned i=0; i<mystash->getNumberOfStoredValues();++i){
        myvalues[i].resize(NumberOfQuantities);
        mystash->retrieveSequentialValue(i, false, myvalues[i] );
    }
    NumberOfDerivatives=amc->getNumberOfDerivatives();
    MultiValue val_k( amc->getNumberOfQuantities(), amc->getNumberOfDerivatives() ); val_k.clearAll();
    myder.resize(NumberOfStoredValues);
    for(unsigned k=0; k<NumberOfStoredValues; ++k) {
        stash->retrieveDerivatives( k, false, val_k );
        myder[k].resize(NumberOfQuantities);
        for(unsigned i=0; i<NumberOfQuantities; ++i){
            myder[k][i].resize(NumberOfDerivatives);
            for(unsigned j=0; j<NumberOfDerivatives; ++j) myder[k][i][j]=val_k.getDerivative(i,j);
        }
    }
}else log<<"*** No ActionWithVessel ***\n";
...
}

Of course, several new cmd_... are defined.

void PlumedMain::cmd(const std::string & word,const TypesafePtr & val) {
...
    case cmd_getNumberOfTasks:
        CHECK_INIT(initialized,word);
        CHECK_NOTNULL(val,word);
        *val.get<int*>()=getNumberOfTasks();
        break;
   case cmd_getNumberOfQuantities:
        CHECK_INIT(initialized,word);
        CHECK_NOTNULL(val,word);
        *val.get<int*>()=getNumberOfQuantities();
        break;
   case cmd_getMultiColvars:
        CHECK_INIT(initialized,word);
        plumed_assert(nw==2);
        CHECK_NOTNULL(val,word);
        int itask; Tools::convertNoexcept(words[1], itask);
        *( val.template get<const double**>())=&myvalues[itask][0];
        break
   case cmd_getDerivativeOfMCV:
        CHECK_INIT(initialized,word);
        plumed_assert(nw==3);
        CHECK_NOTNULL(val,word);
        getDerivativesOfMCV(val, words[1], words[2]);
        break;
...
}

where getDerivativesOfMCV is like this:

void PlumedMain::getDerivativesOfMCV(const TypesafePtr& val, const std::string& word1, const std::string& word2){
  int itask;
  int jtask;
  Tools::convertNoexcept(word1, itask);
  Tools::convertNoexcept(word2, jtask);
  auto vval=val.template get<const double**>();
  if(!myder.empty()){
    *vval=&myder[itask][jtask][0];
  }else *vval=NULL;
}

When calling Plumed from LAMMPS, fix_plumed.cpp to be specific, I changed p->cmd("performCalc") to p->cmd("performCalcNoUpdate").

tuoping avatar Nov 08 '21 09:11 tuoping