Question about CustomCentroidBondForce

[Kailong Mao]

Hi Friends:

Good afternoon! I added a CustomCentroidBondForce between two groups of atoms and apparently it had no effect during the simulation. The two groups of atoms deviated very far away from each other after 500000 steps. I then hard-coded the energy expression to be "10000", added the CustomCentroidBondForce to system, and had OpenMM output the potential energy at step 0. I compared the potential energy calculated with and without this extra force and found out that, apparently, adding this force had no effect on the potential energy.

The code for adding

Code: Select all

OpenMM_CustomCentroidBondForce* force = OpenMM_CustomCentroidBondForce_create(2, "100000");
...
OpenMM_IntArray* bondGroups = OpenMM_IntArray_create(0);
OpenMM_IntArray_append(bondGroups, 0);
OpenMM_IntArray_append(bondGroups, 1);
OpenMM_DoubleArray* bondParameters = OpenMM_DoubleArray_create(0);
OpenMM_CustomCentroidBondForce_addBond(force, bondGroups, bondParameters);

OpenMM_System_addForce(system, (OpenMM_Force *) force);

I tested the code on both the Reference and CUDA platforms and the results were the same.

Thank you very much!

[Peter Eastman]

An energy expression of "100000" means that the energy is always exactly 100000, independent of the positions of the groups. So it's not going to apply any force at all! The force is the negative gradient of the energy, so if you want it to apply forces to the groups, the energy has to depend on their positions. For example, if you set the energy expression to "0.5*k*distance(g1,g2)^2", they'll attract each other with a force proportional to the distance between them.


Peter

[Kailong Mao]

Hi Peter:

Thank you for the quick response! I hard-coded the energy expression because I wanted to compare the potential energy at step 0 with and without the CustomCentroidBondForce. Anyway, I just recompiled my code using the energy expression "100000*distance(g1,g2)^2". I compared the potential energy with and without the CustomCentroidBondForce at step 0. And there is still no difference between the two.

Thanks!

[Peter Eastman]

Could you post your complete code then? We'll need more information to see exactly what you're doing.

Peter

[Kailong Mao]

No problem! Thank you very much! This code actually sits in the OpenMM-Tinker interface file. So it is pretty bulky. But the only lines I added to this file are:

Code: Select all

static OpenMM_IntArray* getGroup(int* kgrp, int* igrp, int idx){
   int start = igrp[2 * idx] - 1;
   int end = igrp[2 * idx + 1] - 1;
   OpenMM_IntArray* group = OpenMM_IntArray_create(0);
   for (int i = start; i < end + 1; i++){
      OpenMM_IntArray_append(group, kgrp[i] - 1);
   }
   printElements(group);
   return group;
}

static OpenMM_DoubleArray* getWeights(int* kgrp, int* igrp, int idx, OpenMM_System* system){
   int start = igrp[2 * idx] - 1;
   int end = igrp[2 * idx + 1] - 1;
   OpenMM_DoubleArray* weights = OpenMM_DoubleArray_create(0);
   for (int i = start; i < end + 1; i++){
      OpenMM_DoubleArray_append(weights, OpenMM_System_getParticleMass(system, kgrp[i] - 1));
   }
   fprintf (stdout, "\n Masses of atoms:");
   printDoubles(weights);
   return weights;
}

static void setupCentroidConstraints(OpenMM_System* system, FILE* log) {
   // In the expression below, u is the upper bound, l is lower bound
   OpenMM_CustomCentroidBondForce* force = OpenMM_CustomCentroidBondForce_create(2, 
                               "step(distance(g1,g2)-u)*k*(distance(g1,g2)-u)^2+step(l-distance(g1,g2))*k*(distance(g1,g2)-l)^2");
   OpenMM_CustomCentroidBondForce_addPerBondParameter(force, "k");
   OpenMM_CustomCentroidBondForce_addPerBondParameter(force, "l");
   OpenMM_CustomCentroidBondForce_addPerBondParameter(force, "u");

   // Add all the groups to force
   for (int i = 1; i < group__.ngrp + 1; i++){
      OpenMM_CustomCentroidBondForce_addGroup(force, 
                                              getGroup(group__.kgrp, group__.igrp, i),
                                              getWeights(group__.kgrp, group__.igrp, i, system));
   }
   fprintf (stdout, "\n #groups: %d \n", OpenMM_CustomCentroidBondForce_getNumGroups(force));


   for (int i = 0; i < restrn__.ngfix; i++) {
      OpenMM_IntArray* bondGroups = OpenMM_IntArray_create(0);
      OpenMM_IntArray_append(bondGroups, restrn__.igfix[2 * i] - 1);
      OpenMM_IntArray_append(bondGroups, restrn__.igfix[2 * i + 1] - 1);
      printElements(bondGroups);

      OpenMM_DoubleArray* bondParameters = OpenMM_DoubleArray_create(0);
      OpenMM_DoubleArray_append(bondParameters, restrn__.gfix[2 * i]);
      OpenMM_DoubleArray_append(bondParameters, restrn__.gfix[2 * i + 1]);
      OpenMM_DoubleArray_append(bondParameters, restrn__.gfix[2 * i + 2]);
      printDoubles(bondParameters);

      OpenMM_CustomCentroidBondForce_addBond(force, bondGroups, bondParameters);
      
   }
   fprintf (stdout, "\n Restraints: Added all constraints\n");
   OpenMM_System_addForce(system, (OpenMM_Force *) force);
}
...
setupCentroidConstraints(omm->system, log);

Please see attached for the complete file. Also, please let me know if you need all the Fortran code to compile this code. I know the atom indices, atom weights (mass), and parameters - k, l, u are read into this part of the code correctly. I confirmed this with all the print statements in the code.

Thank you very much!

[Kailong Mao]

Hi Peter:

It looks like something is wrong with the OpenMM installed on our computer. I wrote a Python script using CustomCentroidBondForce and apparently it works well on another machine. Thank you very much for your help!