csg_stat - missing exclusions in polyaromatic rings (& branched polymers?)

[GoogleCodeExporter]

Sharp peaks appear in RDFs calculated with csg_stat; due to intra-chain pairs 
that should have been excluded from non-bonded RDF calculation. Notably, g_rdf 
(of Gromacs) does the job correctly - no such peaks for the same 
system/topology (see the provided example).

The issue does not occur for linear or two tail polymer chains, e.g. in lipid 
molecules no intra-chain peaks are seen even with nrexcl=3. That is, in simpler 
than above cases changing nrexcl affects the RDFs in a (seemingly) correct way.

What steps will reproduce the problem?

1. Attached is an example of the system where the issue is present for a number 
of non-bonded pairs, most noticeable for B-B CG beads (see the distribution 
files).

2. It appears that changing nrexcl paramater (next to the molecule name in the 
.top file) does not affect the RDF results of csg_stat.

3. It does not matter if csg_stat called for CG or full-atom (without solvent) 
system; the example includes files for both representations, with short .xtc 
trajectories.

What is the expected output? What do you see instead?

- The expected output is given in B-B.dist.g_rdf-no_peaks.xvg as obtained using 
g_rdf

Original issue reported on code.google.com by abruk...@googlemail.com on 10 Oct 2014 at 11:33

Attachments:

• 40TP6-exclusions-issue.zip

[GoogleCodeExporter]

What is the command to get B-B.dist.new-csg_stat-peaks using csg_stat?

Original comment by christop...@gmail.com on 10 Oct 2014 at 11:43

[GoogleCodeExporter]

The standard csg_stat call, say this:

$ csg_stat --nt 4 --top 40TP6-fa-no_sol.tpr --trj 40TP6-fa-no_sol.xtc --options 
40TP6-cg_int_map.xml --cg 40TP6-cg_map.xml

or this (produces the same output):

csg_stat --nt 4 --top 40TP6-cg.tpr --trj 40TP6-cg.xtc --options 
40TP6-cg_int_map.xml --cg 40TP6-cg_map11.xml

Original comment by abruk...@googlemail.com on 10 Oct 2014 at 11:53

[GoogleCodeExporter]

Sorry, forgot to mention: substitute the .xtc file names for the ones provided, 
but beware - the latter are very short! the original .xtc files cannot be 
supplied due to their huge size (obviously).

Note that the command(s) will produce ALL distributions for the system, 
including B-B.dist.new where you will see the sharp peaks (peaks will show up 
in a few distributions for the beads in polyaromatic core/rings).

Original comment by abruk...@googlemail.com on 10 Oct 2014 at 11:58

[GoogleCodeExporter]

...use vmd for visualization of the system! :-)

Original comment by abruk...@googlemail.com on 10 Oct 2014 at 12:00

[GoogleCodeExporter]

Ok, I think I got it, this is a feature not a bug.

With "--cg" a mapping is performed, so whatever the exclusions on the atomistic 
level were does not matter as there is no generic way to map exclusions. (Also 
one wants to have control over the exclusions on the cg level.)

Exclusions after a mapping are set from the bonded interactions defined in the 
mapping file. Convention: all bonds, angles and dihedrals are excluded. 
$ csg_dump --top 40TP6-cg.tpr --cg 40TP6-cg_map11.xml --excl
will show the exclusions VOTCA generates.
To add extra exclusions, you will need add some "phantom" bonds.

Also see the following link for details:
<https://groups.google.com/d/topic/votca/NDvv-iPsL0I/discussion>

A nrexcl option in file would be nice - patches welcome!

Original comment by christop...@gmail.com on 10 Oct 2014 at 12:15

• Changed state: WontFix

[GoogleCodeExporter]

Christoph,

What you say implies nrexcl in CG topology is not recognised (not taken into 
account) by VOTCA tools at all. It is so?

I just wonder why g_rdf does the job correctly (nrexcl=3 which is the standard 
way of dealing with pairs within bonds, angles & dihedrals as you described).

Original comment by abruk...@googlemail.com on 10 Oct 2014 at 12:20

[GoogleCodeExporter]

Yes, if you map ("--cg" option) nrexcl from tpr is ignored by design. Without 
mapping nrexcl from tpr is used.

nrexcl=3 is the standard for some molecules, but surely not for all and hence 
we didn't want to hard-code it in VOTCA. In some sense VOTCA is using nrexcl=1. 
nrexcl!=1 can easily be achieved by adding extra "phantom" bonds/angles.

Original comment by christop...@gmail.com on 10 Oct 2014 at 12:29

[GoogleCodeExporter]

OK, I see now the root of the problem - csg_stat does not have --no-map option!
Otherwise I would use that, of course.

Thanks for the tip on phantom bonds. Will try.

Original comment by abruk...@googlemail.com on 10 Oct 2014 at 12:49

[GoogleCodeExporter]

Hi,

I am not sure what do you mean by '--no-map'? Do you mean you do not want to 
map and then determine the RDFs of mapped CG beads? csg_stat does allow to 
determine rdfs from original (non-mapped) trajectory, for that you just run 
csg_stat with --cg option. With no '--cg' option, as Chirstoph pointed out, 
csg_stat will use exclusion ruled from the reference .tpr file.

---
Sikandar

Original comment by masha...@votca.org on 13 Oct 2014 at 7:34

[GoogleCodeExporter]

Hi Sikandar,

First, thanks to Christoph's hints (partly by email) I have by now figured out 
how to obtain the correct RDFs by using nrexcl in Gromacs .top ( and hence 
.tpr) for the beforehand CG-mapped system. 

Second, csg_map has option --no-map, but csg_stat only has --cg, but can be run 
without this option which appears to be equivalent --no-map (instead of it). It 
is of course confusing, but as soon as one learns this, it's alright.

To resume, this issue is probably about making the use of csg_stat consistent 
with csg_map (in terms of options), on the one hand, but also there is an extra 
confusion:

As of now, in order to obtain correct RDFs (using correct nrexcl) csg_stat has 
to be used without mapping (skipping --cg option and running it for the 
CG-mapped trajectory), but to obtain bonded distributions one has to always use 
--cg option (with fictitious 1:1 mapping, if using CG-mapped traj), because 
otherwise csg_stat cannot get any info about the topology of bonding.

So I guess, this is actually a low priority fix, but it would make using 
csg_stat more transparent and consistent.

Cheers / Andrey

Original comment by abruk...@googlemail.com on 14 Oct 2014 at 1:55

[GoogleCodeExporter]

I made the nrexcl option a feature request (issue #163).

The name csg_map implies that is does a mapping by default, while csg_stat 
doesn't imply that. Hence csg_map as an option to disable the default mapping 
and csg_stat not.

Original comment by christop...@gmail.com on 14 Oct 2014 at 3:09