bug in communication of mass between NAMD2.9 and PLUMED2.1

[Glen Hocky]

Hi Everyone, 

I've discovered a small but dangerous bug when using NAMD with plumed (I have only rigorously tested with NAMD 2.9 and the latest development version of PLUMED). It seems that the masses are not being passed correctly from NAMD to PLUMED, resulting in incorrect calculations when using the COM (or CENTER ... MASS) keyword. Having done some digging, I have found that the masses array is initially filled with the correct masses for that atoms defined in groups, but then the masses array is later filled in with sequential masses for the atoms in the topology of the system (i.e. for atoms 1-n, with n the number of requested atoms from all of the Actions).

Interestingly, this problem can be fixed by accessing all of the atoms needed in a sequential manner just a single time, e.g. I have found that doing something like this

DUMPATOMS STRIDE=5000000 FILE=out.xyz ATOMS=1-5000

will fix the problem where the maximum atom index needed is < 5000

Although there is a way around the problem, it would be great if we could get this fixed. Information on my debugging efforts are below.

Best,

Glen

--

I have attached a demonstration of this problem, where I have taken the inputs for the NAMD tutorial simulation of ubiquitin in a water sphere, and call a broken and "working" plumed file, which computes the COM of the alpha carbon and nitrogen for residues 1 and 3 as the first group, and 2 and 4 as the second group, and then computes the distance between these.

Both simulations run 0 time steps, to demonstrate that the problem occurs even before dynamics

To run:

The following commands will produce xyz files showing the correct and incorrect COM positions, and a correct and incorrect distance calculation

namd2 ubq_ws_eq_plumed-broken.conf

namd2 ubq_ws_eq_plumed-works.conf

I have also included an independent check, a NAMD tcl script which computes the COMs and distances, and agrees with the "working" plumed calculation

To run:

vmd -e measure_cv.tcl -dispdev none

Please let me know if you have any problems here

--

More about debugging. I have done some manual debugging via printing out various variables during the example simulations above. The problem seems to be related to the plumed-patch share() function in  NAMD_2.9_Source/src/ComputeMgr.C.

In the example simulations that I have attached, this function is called 2x. The first time, it calls createFullList and then fills the masses array correctly. Some of my print outs are here

in plumed createFullList - 0

in plumed createFullList - 4

in plumed createFullList - 19

in plumed createFullList - 21

in plumed createFullList - 36

in plumed createFullList - 38

in plumed createFullList - 55

in plumed createFullList - 57

full list created - n: 8

get full list run: 0 - n: 8

in redo - n: 8

0:14.007000 4:12.011000 19:14.007000 21:12.011000 36:14.007000 38:12.011000 55:14.007000 57:12.011000 

Before share() is called again, I see the first computation of the COMs, and already the masses accessed by COM are not correct

TCL: Running for 0 steps

Computing COM - natoms: 4

mass 0: 14.007000

mass 1: 1.008000

mass 2: 12.011000

mass 3: 1.008000

Computing COM - natoms: 4

mass 0: 1.008000

mass 1: 1.008000

mass 2: 12.011000

mass 3: 1.008000

Finally, when share() is called again I see

full list created - n: 0

get full list run: 8 - n: 0

// PRINT OUT full masses array

i 0 mass 14.007000 

i 1 mass 1.008000 

i 2 mass 1.008000 

i 3 mass 1.008000 

i 4 mass 12.011000 

i 5 mass 1.008000 

i 6 mass 12.011000 

i 7 mass 1.008000 

in redo - n: 0

///

However, in the case where I use the DUMPATOMS command to access atoms 1-100

At this final step, the masses array has 100 elements with masses of atoms 1-100. Presumably this is related to why the COM command is able to access the correct masses for the 8 selected atoms in the groups

[Giovanni Bussi]

Hi

many thanks for pointing this out!!

I am not familiar at all with namd, but I had a look at the patch and it could be that the problem is related to these lines:

+    for(int i=0;i<index.size();i++){

+      Vector coord;

+      getPosition(index[i],coord);

+      positions[3*i+0]=coord.x;

+      positions[3*i+1]=coord.y;

+      positions[3*i+2]=coord.z;

+      masses[i]=Node::Object()->molecule->atommass(i);

+      charges[i]=Node::Object()->molecule->atommass(i);

+    };

Can you just try to change to

+    for(int i=0;i<index.size();i++){

+      Vector coord;

+      getPosition(index[i],coord);

+      positions[3*i+0]=coord.x;

+      positions[3*i+1]=coord.y;

+      positions[3*i+2]=coord.z;

+      masses[i]=Node::Object()->molecule->atommass(index[i]);

+      charges[i]=Node::Object()->molecule->atommass(index[i]);

+    };

(change in bold) and check if this fixes the bug?

Notice that, if this is the problem, also charges are wrong.

Thanks a lot!

Giovanni

--
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 post to this group, send email to plumed...@googlegroups.com.
Visit this group at http://groups.google.com/group/plumed-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/07f2a653-af03-42f6-905a-04c66da0679e%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Giovanni Bussi]

By the way, plumed only cares about the masses received at step 1 (they are assumed to be constant). That's the reason why your (smart!) fix works. And notice that since v2.2 they can be printed with DUMPMASSCHARGE on a file so as to double check (or to read them later from the "plumed driver")

Giovanni

[Glen Hocky]

Hi Giovanni,

Thanks for your quick response! I didn't know when/why easy_calc was being called, so i didn't think to check there. Your fix seems to solve the problem!

By the way, for the charge, should it be atomcharge as below?

+      masses[i]=Node::Object()->molecule->atommass(index[i]);

+      charges[i]=Node::Object()->molecule->atomcharge(index[i]);

You received this message because you are subscribed to a topic in the Google Groups "PLUMED users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/plumed-users/rsBD9JoOtN4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to plumed-users...@googlegroups.com.

To post to this group, send email to plumed...@googlegroups.com.
Visit this group at http://groups.google.com/group/plumed-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAPLm8ZKs627NAELBXkApPy1nDKu7ghBEX0xK5n94gKn87B6jfw%40mail.gmail.com.

[Giovanni Bussi]

Arg, you are obviously right!

I opened a github issue on this github.com/plumed/plumed2/issues/162 with some comment about why the bug was not noticed before. I think the error has always been there but the bug was only visible at first step in v2.0. Since v2.1 the bug affects the entire simulation.

I don't want to install namd, so let me know if the fix is fine I will add it as is to the repository (luckily it can be applied without redoing the diff).

Thanks again for the bug and the solution!

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAGMfakkqrabPBEBj-B0Gett-nk6J8OQtuaaO3eHEi4eN0TZ68w%40mail.gmail.com.

[Glen Hocky]

As far as my example can detect, this fix solves the problem for the masses.

I tried to check the charges using the DIPOLE command, but I'm getting the result: "charges were not passed to plumed". If charge CV's aren't enabled for NAMD, perhaps that's why no one noticed the charge bug before.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAPLm8Z%2BN7tO5iT%3D82UiRNN5%3Dj88mwpTHSM9OcEXaMUyT8maXAQ%40mail.gmail.com.

[Giovanni Bussi]

Perhaps you should also change:

+    cmd("setMasses",&masses[0]);

+    cmd("setPositions",&positions[0]);

+    cmd("setForces",&forces[0]);

to

+    cmd("setMasses",&masses[0]);

+    cmd("setCharges",&charges[0]);

+    cmd("setPositions",&positions[0]);

+    cmd("setForces",&forces[0]);

Thanks!

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAGMfakkr2VTiy%2BL0Tjq7qS-Lmtzhww7m4r1A76fLHVtt2ZFYOw%40mail.gmail.com.

[Glen Hocky]

Yes, now the charges are set and seem correct!

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAPLm8ZLA0bvJOX4g8y2X5GEDBbiUUwE1hVYVySZqvk7P_kKT9Q%40mail.gmail.com.