Re: [plumed-users] Trouble testing new CV - Julich Tutorial

[DAVID SILVA SANCHEZ]

Dear Gareth,

First, I am really sorry for the confusion with the e-mails. I erased the message you answered to, because I wanted to write it more clearly, but now I can't answer you directly, So I just created a new thread.

Second, thank you for your help, that fixed it!

However, I still have troubles with the results I obtain with my code, basically my CV is 0 for every frame. Is there some kind of solution to this tutorial? Here is my compute function (I also attached my cpp and plumed files):

double SecondShell::compute(const unsigned &tindex, AtomValuePack &myatoms) const

{

// Calculate the coordination number

double dfunc, hb, d;

for (unsigned i = 1; i < myatoms.getNumberOfAtoms(); ++i)

{

Vector &distance = myatoms.getPosition(i);

double d2;

if ((d2 = distance[0] * distance[0]) < rcut2 &&

(d2 += distance[1] * distance[1]) < rcut2 &&

(d2 += distance[2] * distance[2]) < rcut2 &&

d2 > epsilon)

{

d = std::sqrt(d2);

//Calculate the value of the integral

hb = histogramBead.calculate(d, dfunc);

//acummulate it

accumulateSymmetryFunction(1, i, hb, (dfunc)*distance, (-dfunc) * Tensor(distance, distance), myatoms);

}

}

return myatoms.getValue(1);

}

I am not really sure what the problem is. I tried multiplying the result of the integral by sqrt(2*pi), since in the tutorial the integral is not normalized, but that didn't fix it.

Note: I checked using different boundaries and the CV is not zero when the boundaries are big (e.g, 1-59 instead of 57-59).

[Gareth Tribello]

Hi David

Here is how I did it:

class SecondSphere : public MultiColvar {

private:

  double rcut2;

  HistogramBead bead;

public:

  static void registerKeywords( Keywords& keys );

  explicit SecondSphere(const ActionOptions&);

// active methods:

  virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;

/// Returns the number of coordinates of the field

  bool isPeriodic(){ return false; }

};

PLUMED_REGISTER_ACTION(SecondSphere,"SECONDSPHERE")

void SecondSphere::registerKeywords( Keywords& keys ){

  MultiColvar::registerKeywords( keys );

  keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB");

  // Use actionWithDistributionKeywords

  keys.add("compulsory","BEAD","the histogram we are using");

  keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("MAX");

  keys.use("MIN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");

  keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");

}

SecondSphere::SecondSphere(const ActionOptions&ao):

PLUMED_MULTICOLVAR_INIT(ao)

{

  std::string mybeadstr, errors; parse("BEAD",mybeadstr);

  bead.set( mybeadstr, errors ); bead.isNotPeriodic();

  if( errors.size()!=0 ) error( errors );

  // Read in the atoms

  int natoms=2; readAtoms( natoms );

  // And setup the ActionWithVessel

  checkRead();

}

double SecondSphere::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {

   double value=0, dfunc;

   // Calculate the coordination number

   double d2, sw;

   for(unsigned i=1;i<myatoms.getNumberOfAtoms();++i){

      Vector& distance=myatoms.getPosition(i);  // getSeparation( myatoms.getPosition(0), myatoms.getPosition(i) );

      double d = distance.modulo();

      sw = bead.calculate( d, dfunc ); dfunc /= d;

      value += sw;

      addAtomDerivatives( 1, 0, (-dfunc)*distance, myatoms );

      addAtomDerivatives( 1, i,  (dfunc)*distance, myatoms );

      myatoms.addBoxDerivatives( 1, (-dfunc)*Tensor(distance,distance) );

   }

   return value;

}

This was multiple years ago, however, and the accumulate symmetry function that you have used looks like it makes sense.

Have you tried printing out the d and hb values that you are calculating and seeing if any of them are non-zero?

Gareth

"La información aquí contenida es para uso exclusivo de la persona o entidad de destino. Está estrictamente prohibida su utilización, copia, descarga, distribución, modificación y/o reproducción total o parcial, sin el permiso expreso de Universidad de Antioquia, pues su contenido puede ser de carácter confidencial y/o contener material privilegiado. Si usted recibió esta información por error, por favor contacte en forma inmediata a quien la envió y borre este material de su computador. Universidad de Antioquia no es responsable por la información contenida en esta comunicación, el directo responsable es quien la firma o el autor de la misma."
UdeA

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/469c3946-16eb-4fa5-bfba-3261fe81e645n%40googlegroups.com.

[DAVID SILVA SANCHEZ]

Dear Gareth,

Thank you for sharing your code with me.

• "Have you tried printing out the d and hb values that you are calculating and seeing if any of them are non-zero?"

Yes, and they are non-zero. Maybe the trajectory and structure used in the regtests changed? Because I don't think my plumed.dat has errors (I hope).

I tried your code, adapted it a bit because it didn't work with the current version (just Multicolvar -> MulticolvarBase) and the same happened. So I think I am going to leave it at that and be happy knowing that I probably programmed it correctly and learned a lot.

I am actually working on a CV that will probably be implemented in plumed, so this discussion is and will be really helpful.

Thank you again,

David