com/chunk and vcm/chunk problems

Dear Lammps users,

I’ve encountered a strange problem when using com/chunk and vcm/chunk computes: defining chunks to be polymer atoms of 3 different polymer chains in water works well for calculation of MSD and gyration radii of corresponding chains/chunks, however center-of-mass and its velocity are always inf-numbers or unphysical values of 10^8. I suspect a bug in the output of these 2 quantities, since MSD, which tracks the value of COM of each chunk looks pretty normal even for long simulation runs over 1 ns. Below is the portion of input script used for test case:

units real

atom_style full

boundary p p p

variable Na equal 1000000

variable T equal 320.0

variable tstep equal 1.0

pair_style lj/long/coul/long long long 10.0

bond_style harmonic

angle_style harmonic

dihedral_style opls

improper_style class2

read_data 3decamer.data

include pairparameters.dat

group poly type 1 2 3 4 5 6 7

group water type 8 9

compute NIP poly chunk/atom molecule nchunk once ids once

---------- Run Minimization ---------------------#

reset_timestep 0

fix 1 all box/relax iso 1.0 vmax 0.001

thermo 100

thermo_style custom step pe etotal lx ly lz press pxx pyy pzz

min_style cg

minimize 1e-6 1e-6 5000 10000

unfix 1

write_restart restartfile.minimized

velocity all create ${T} 4928459 mom yes rot yes dist gaussian

compute RG poly gyration/chunk NIP

compute MSD poly msd/chunk NIP

compute COM poly com/chunk NIP

compute VCM poly vcm/chunk NIP

fix 2 all nvt temp {T} {T} 50

fix PRINT poly ave/time 10 1 10 c_RG c_COM c_VCM c_MSD file polymer.out mode vector

run 1000

the output of the ave/time script is as follows:

0 3

1 6.7764 0 0 0 -inf -inf -inf 0 0 0 0

2 6.58678 0 0 0 inf inf -inf 0 0 0 0

3 7.08112 0 0 0 -inf inf inf 0 0 0 0

10 3

1 6.82073 0 0 0 inf -inf -inf 0.00544238 0.0675754 8.38274e-09 0.0730178

2 6.81958 0 0 0 inf inf inf 0.192354 0.325749 0.015632 0.533736

3 6.80645 0 0 0 -inf inf inf 0.00193617 0.26595 0.0158504 0.283736

20 3

1 6.80703 0 0 0 inf -inf -inf 0.00518961 0.0665801 5.39325e-05 0.0718236

2 6.79353 0 0 0 inf inf inf 0.19582 0.320532 0.0162308 0.532583

3 6.7816 0 0 0 inf inf inf 0.00183996 0.267282 0.0171962 0.286318

30 3

1 6.77937 0 0 0 inf -inf -inf 0.0046755 0.0653146 0.000196425 0.0701866

2 6.75628 0 0 0 inf inf inf 0.199732 0.315468 0.0168438 0.532045

3 6.74668 0 0 0 inf inf inf 0.00143704 0.268844 0.0186708 0.288952

I have tried different initial configurations and restarts from previous runs, still no luck with COM and VCM portion, though MSD and gyration radii look reasonable to me. Any ideas, thoughts?
If needed, the initial configuration can be provided as well.
Best,
Vitalie

Dear Lammps users,

I've encountered a strange problem when using com/chunk and vcm/chunk
computes: defining chunks to be polymer atoms of 3 different polymer chains
in water works well for calculation of MSD and gyration radii of
corresponding chains/chunks, however center-of-mass and its velocity are
always inf-numbers or unphysical values of 10^8. I suspect a bug in the
output of these 2 quantities, since MSD, which tracks the value of COM of
each chunk looks pretty normal even for long simulation runs over 1 ns.
Below is the portion of input script used for test case:

units real

atom_style full

boundary p p p

variable Na equal 1000000

variable T equal 320.0

variable tstep equal 1.0

pair_style lj/long/coul/long long long 10.0

bond_style harmonic

angle_style harmonic

dihedral_style opls

improper_style class2

read_data 3decamer.data

include pairparameters.dat

group poly type 1 2 3 4 5 6 7

group water type 8 9

compute NIP poly chunk/atom molecule nchunk once ids once

please try removing "ids once" and see if that helps in your case.
there were some problems with that flag before, but it looks as if
those are not fully resolved.

axel.

...and the following change to the source code should resolve the
problem when using "ids once".
diff --git a/src/compute_com_chunk.cpp b/src/compute_com_chunk.cpp
index 9219455..0af54d9 100644
--- a/src/compute_com_chunk.cpp
+++ b/src/compute_com_chunk.cpp
@@ -88,8 +88,8 @@ void ComputeCOMChunk::setup()
   // done in setup, so that ComputeChunkAtom::setup() is already called

   if (firstflag && cchunk->idsflag == ONCE) {
- firstflag = massneed = 0;
     compute_array();
+ firstflag = massneed = 0;
   }
}

diff --git a/src/compute_vcm_chunk.cpp b/src/compute_vcm_chunk.cpp
index 0fcb22c..dd2c075 100644
--- a/src/compute_vcm_chunk.cpp
+++ b/src/compute_vcm_chunk.cpp
@@ -88,8 +88,8 @@ void ComputeVCMChunk::setup()
   // done in setup, so that ComputeChunkAtom::setup() is already called

   if (firstflag && cchunk->idsflag == ONCE) {
- firstflag = massneed = 0;
     compute_array();
+ firstflag = massneed = 0;
   }
}

Axel,

your patch solved the problem, as usually. Thanks a lot!
Hope it will be pushed to the next lammps patch.

vitalie

good catch - will be in the next patch.

Steve