Reading Lammps topology and dump into csg_stat

[Daniel Allen]

Hi there,

I am trying to read in a LAMMPS dump file for a simulation of liquid butane into csg_stat as a simple example of how to map and process LAMMPS trajectories, before moving onto my more complex system of long chain polymers.

I know that there are a couple of previous questions regarding LAMMPS topologies but I cannot quite get it working.

First I tried using a .pdb file specified in my topol.xml file but this gave an error:

"an error occurred:
P-P.pdb: unknown topology format" which I find strange because in the manual it gives an example of using a pdb file for topology, any ideas why this error would occur?

Anyway, I then decided to try and use the lammps .dump file as the topology so topol.xml now looks like this:

<topology base="atoms_nvt_prod.dump">/ITEM
  <molecules>
     <clear/>
     <define name="B" first="1" nbeads="14" nmols="100"/>
  </molecules>
</topology>

So I have defined one molecule type which consists of 14 atoms and there are 100 molecules in the system beginning at atom 1.

I think the problem lies in the coarse-grained mapping file. When listing atoms which belong to a particular bead, if the syntax is RESID:RESNAME:ATOMNAME then I'm struggling to work out what they should be for my system.

Here is the top of my dump file:

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
13916
ITEM: BOX BOUNDS pp pp pp
-27.8209 27.7639
-27.7926 27.8726
-27.6615 27.7715
ITEM: ATOMS id type x y z
12765 2 -21.1019 -18.5716 -26.1089
12766 1 -21.1227 -18.563 -24.9982
12768 1 -20.1424 -19.0188 -26.4462
9768 1 -25.3608 -20.7513 -21.6347
9766 2 -25.8077 -20.8022 -20.6188
9762 1 -23.0727 -20.4904 -19.9854
9764 1 -25.494 -18.664 -19.9916
9763 2 -25.1286 -19.6941 -19.7923
9767 1 -26.8648 -20.4951 -20.7686
9769 2 -25.9722 -22.1004 -19.9808
9760 1 -23.565 -19.281 -21.2077
9759 2 -23.6425 -19.548 -20.132
13105 2 -25.5908 -19.0274 -24.4098
638 2 -22.4319 -26.8526 -26.9045
9772 1 -26.1693 -22.9545 -20.6634
13106 1 -25.5393 -19.9713 -24.9935
11466 1 -22.2057 -26.6297 -23.7546
11454 1 -26.0494 -24.6427 -23.5225
11458 1 -24.8424 -26.6034 -24.4689
3362 1 -20.4076 -20.3751 -20.0356
639 1 -22.2076 -27.6125 -26.1257

etc ...

Would RESID just be 1? There is no information about residues in the dump file, does there need to be? It is possible to printout the molecule number in the dump file.

For RESNAME, I wouldn't know what to use either.

For ATOMNAME I could use the atom type, so atom name would be "2" for C and "1" for H?

The choice for these labels would be more obvious from my .pdb file:


HETATM    1  C7   B   A   1      -5.985   7.723  -2.9591.00  0.00           C
HETATM    2  H71  B   A   1      -5.479   7.242  -2.1481.00  0.00           H
HETATM    3  H72  B   A   1      -5.406   7.617  -3.8531.00  0.00           H
HETATM    4  C8   B   A   1      -6.944   7.270  -3.1001.00  0.00           C
HETATM    5  C9   B   A   1      -7.523   7.376  -2.2061.00  0.00           C
HETATM    6  H81  B   A   1      -6.818   6.231  -3.3221.00  0.00           H
HETATM    7  H82  B   A   1      -7.450   7.750  -3.9111.00  0.00           H
HETATM    8  H91  B   A   1      -7.017   6.896  -1.3951.00  0.00           H
HETATM    9  C10  B   A   1      -8.482   6.923  -2.3471.00  0.00           C
HETATM   10  H92  B   A   1      -7.648   8.415  -1.9831.00  0.00           H
HETATM   11  H101 B   A   1      -8.357   5.884  -2.5691.00  0.00           H
HETATM   12  H102 B   A   1      -8.988   7.404  -3.1581.00  0.00           H
HETATM   13  H73  B   A   1      -6.112   8.781  -2.7321.00  0.00           H
HETATM   14  H103 B   A   1      -9.071   7.031  -1.4361.00  0.00           H

So for the first atom it would be   1:B:C7 for example?

 

Any help would be greatly appreciated,

Thanks, Dan.

[Christoph Junghans]

Hi Dan,

welcome.

2014-10-20 4:37 GMT-06:00 Daniel Allen <daniela...@gmail.com>:
> Hi there,
>
> I am trying to read in a LAMMPS dump file for a simulation of liquid butane
> into csg_stat as a simple example of how to map and process LAMMPS
> trajectories, before moving onto my more complex system of long chain
> polymers.
>
> I know that there are a couple of previous questions regarding LAMMPS
> topologies but I cannot quite get it working.
>
> First I tried using a .pdb file specified in my topol.xml file but this gave
> an error:
>
> "an error occurred:
> P-P.pdb: unknown topology format" which I find strange because in the manual
> it gives an example of using a pdb file for topology, any ideas why this
> error would occur?

Which version of VOTCA is that?
In 1.2 the pdb is based on the gromacs reader, which means without
Gromacs you won't be able to use pdb files.
However in 1.3 the pdb reader is independent from the Gromacs library,
so make sure you have the 1.3-dev version installed before going on:
<https://code.google.com/p/votca/wiki/Installing#Development_Version>
Also the lammps dump reader in 1.3 is a bit more versatile.

You can simply run
$ csg_dump --top topol.xml
to see how VOTCA see and how it names the beads.
The residue name/number doesn't really matter as long as the triple
"RESID:RESNAME:ATOMNAME" is unique for every bead.

Cheers,

Christoph

>
>
> Any help would be greatly appreciated,
>
> Thanks, Dan.
>
>
>
>
>

> --
> You received this message because you are subscribed to the Google Groups
> "votca" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to votca+un...@googlegroups.com.
> To post to this group, send email to vo...@googlegroups.com.
> Visit this group at http://groups.google.com/group/votca.
> For more options, visit https://groups.google.com/d/optout.



--
Christoph Junghans
Web: http://www.compphys.de

[Daniel Allen]

Hi Christoph,

thanks for your reply.

So I installed the 1.3-dev version of votca after your response. However, running the command:

$ csg_dump --top topol.xml

resulted in the output:
"an error occurred:
bad lexical cast: source type value could not be interpreted as target"

Have you seen this error before, or can you think of anything obvious which might be causing it?

Cheers, Dan.

[Christoph Junghans]

2014-10-31 5:58 GMT-06:00 Daniel Allen <daniela...@gmail.com>:
> Hi Christoph,
>
> thanks for your reply.
>
> So I installed the 1.3-dev version of votca after your response. However,
> running the command:
>
> $ csg_dump --top topol.xml
>
> resulted in the output:
> "an error occurred:
> bad lexical cast: source type value could not be interpreted as target"
>
> Have you seen this error before, or can you think of anything obvious which
> might be causing it?

This is a boost parsing error. Can you post your topol.xml?

[Daniel Allen]

Sure, topol.xml looks like this:

<topology base= "B-B.pdb">

  <molecules>
     <clear/>
     <define name="B" first="1" nbeads="14" nmols="100"/>
  </molecules>
</topology>

Cheers, Dan.

[Christoph Junghans]

2014-10-31 9:57 GMT-06:00 Daniel Allen <daniela...@gmail.com>:
> Sure, topol.xml looks like this:
>
> <topology base= "B-B.pdb">
> <molecules>
> <clear/>
> <define name="B" first="1" nbeads="14" nmols="100"/>
> </molecules>
> </topology>

Looks ok to me, can you send topol.xml and B-B.pdb to de...@votca.org
and I will have a further look.

Christoph

[Christoph Junghans]

2014-10-31 10:56 GMT-06:00 Christoph Junghans <jung...@votca.org>:
> 2014-10-31 9:57 GMT-06:00 Daniel Allen <daniela...@gmail.com>:
>> Sure, topol.xml looks like this:
>>
>> <topology base= "B-B.pdb">
>> <molecules>
>> <clear/>
>> <define name="B" first="1" nbeads="14" nmols="100"/>
>> </molecules>
>> </topology>
> Looks ok to me, can you send topol.xml and B-B.pdb to de...@votca.org
> and I will have a further look.

Thanks, I got it.

Your pdb is mis-formated. In a pdb file the atom name is supposed to
be column 13-16
(http://deposit.rcsb.org/adit/docs/pdb_atom_format.html#ATOM),
but in your file it is 14-17. Removing the space in column 13 using:
$ sed -i 's/\(.\{12\}\)./\1/' B-B.pdb
fixes that.

However, there was also a try and catch block missing in our code,
which I have fixed now:
https://code.google.com/p/votca/source/detail?r=31117aa998ced63f88ae91dab44146bb93aefbce&repo=csg

Christoph

[Vitalie Botan]

Hi Christoph,

for running IBI with Lammps (or any other supported sim. program) there is a special tag <topol> in the settings.xml for csg_inverse. As I got from the description, this tag passes the topology file to the csg_stat script. However, csg_stat accepts also mapping and bonding information via --cg option. Is there any tag to include mapping/bonding option for csg_stat in the settings.xml of csg_inverse? Or how otherwise excluded interactions will be taken into account?

Best,

--Vitalie

[Christoph Junghans]

Hi,

Yes, the tag is called cg.inverse.map! See the hexane/ibi_bonded
example, where hexane_cg.xml is a simple 1:1 map to add extra bonded
interaction and hence exclusions.

Christoph


>
> Best,
> --Vitalie

[Vitalie Botan]

On Tuesday, May 19, 2015 at 4:10:07 PM UTC+2, Christoph Junghans wrote:

Hi,

Yes, the tag is called cg.inverse.map! See the hexane/ibi_bonded
example, where hexane_cg.xml is a simple 1:1 map to add extra bonded
interaction and hence exclusions.

Christoph

Thanks!
--Vitalie

[Vitalie Botan]

It seems really weird, but this tag is somehow ignored by the csg_inverse script:

<inverse>

   <!-- 260*0.0019872041 lammps units -->

   <kBT>0.516673</kBT>

   <!-- use lammps as simulation program -->

   <program>lammps</program>

   <!-- lammps specific options -->

   <lammps>

     <!-- lammps script to run !-->

     <script>in.lmp</script>

     <!-- topology to be used by  csg_stat !-->

     <topol>nipam_cg.xml</topol>

     <!-- traj file created by lammps !-->

     <traj>traj.dump</traj>

   </lammps>

   <map>map.xml</map>

   <initial_configuration>maindir</initial_configuration>

 

So mapping file is included, but I still get this error:

 ERROR:                                                                                                                                                                                          

# critical: 'csg_stat --nt 4 --options /Volumes/Home/Work/Molstruct/data/decamer/3decamer/T260K/coarse-graining/settings.xml --top nipam_cg.xml --trj traj.dump --begin 0 --first-frame 0' failed #

Apparently, csg_stat is called without --cg map.xml , which breaks the type definitions and all the rest of the mapping.

Any ideas?

--Vitalie

[Christoph Junghans]

You might have to add a dummy bonded interaction in your settings xml
file to enable that feature.

Christoph

[Vitalie Botan]

The problem is that the name of the bonded interaction has to match the name defined in mapping file. Is there any way to disable all calls to csg_stat and other scripts in the case of dummy bonded interactions?

--Vitalie

[Christoph Junghans]

I never though about using a mapping file without bonded action, but
your case actually makes sense!

I changed that in the code:
<https://code.google.com/p/votca/source/detail?r=f2366edff0c74460899c81ade4ee247ee1353e54&repo=csg>

[Vitalie Botan]

It works now! Thanks a lot, that is really immediate response and support ever seen ))

--Vitalie

[Vitalie Botan]

The results, however, look really odd. In my mapping file A-A beads are always 1-3 neighbors apart from each, while A-B and B-B are bonded directly. Those delta-like peaks in the RDFs could possibly result from not being excluded properly?

[Christoph Junghans]

2015-05-21 8:35 GMT-06:00 Vitalie Botan <vitali...@gmail.com>:

The results, however, look really odd. In my mapping file A-A beads are always 1-3 neighbors apart from each, while A-B and B-B are bonded directly. Those delta-like peaks in the RDFs could possibly result from not being excluded properly?

It seems like there maybe some  exclusions missing. Remember VOTCA used Gromacs-like exclusion scheme with nrexcl=1, meaning all atom involved in an bond, angle or dihedral are excluded.

You can check the exclusions with 

$ csg_dump --top XXX --cg YYY --excl

and add more bonds to the mapping file to exclude more atom.

Christoph

--
You received this message because you are subscribed to the Google Groups "votca" group.
To unsubscribe from this group and stop receiving emails from it, send an email to votca+un...@googlegroups.com.
To post to this group, send email to vo...@googlegroups.com.
Visit this group at http://groups.google.com/group/votca.
For more options, visit https://groups.google.com/d/optout.

[Vitalie Botan]

I have checked the exclusions and they seem to be ok, at least it works pretty well with the same mapping on the atomistic trajectory. I suspect it's purely numerical stability and convergence problems. In that respect I would like to ask you how the convergence criteria is actually defined in VOTCA, which functional form it takes (I couldn't find this information in the manual) and what is that damping parameter  (eq. 2.15 of the manual) when doing updates, can it be controlled in the same spirit as Picard iteration scheme? I'm not sure what is the best way to iterate over multicomponent system: currently I update only single potential every single step, so that means that for one iteration only one convergence parameters should be checked. Is there also a way to control iterations with this parameter, i.e. to reject those updates, which lead to bigger convergence parameter? 

--Vitalie

[Christoph Junghans]

2015-05-29 12:57 GMT-06:00 Vitalie Botan <vitali...@gmail.com>:
> I have checked the exclusions and they seem to be ok, at least it works
> pretty well with the same mapping on the atomistic trajectory. I suspect
> it's purely numerical stability and convergence problems. In that respect I
> would like to ask you how the convergence criteria is actually defined in
> VOTCA, which functional form it takes (I couldn't find this information in
> the manual) and what is that damping parameter (eq. 2.15 of the manual)
> when doing updates, can it be controlled in the same spirit as Picard
> iteration scheme? I'm not sure what is the best way to iterate over
> multicomponent system: currently I update only single potential every single
> step, so that means that for one iteration only one convergence parameters
> should be checked. Is there also a way to control iterations with this
> parameter, i.e. to reject those updates, which lead to bigger convergence
> parameter?

Eq 2.15 refers to scaling the update with a constant factor. This
factor and also the update order for multicomponent system is very
specific to the system. I would either update all interactions at the
same time and scale them with 1/n or update one at the time and don't
scale them at all.

We have recently extended the post-add convergence calculation (see
postadd_convergence.sh for details), one can now calculated the
weighted sum of multiple difference of e.g. distribution, potential
with respect to the target or the last staep using L1 or L2 norm.
(Defaults are L1 norm of dist.new - dist.tgt with weight 1).
Using the post-add accumulated convergence script one can track the
convergence over steps. See settings.xml of the urea-water/kb-ibi
tutorial for details how to use that.

To do what you proposed above you will need to write you own custom
post-add script (see the customization section of the manual), but the
overwrite post-add script (postadd_overwrite.sh) seems like a good
starting point.

Christoph

>
> --Vitalie
>
> On Thursday, May 21, 2015 at 7:45:33 PM UTC+2, Christoph Junghans wrote:
>>
>>
>> It seems like there maybe some exclusions missing. Remember VOTCA used
>> Gromacs-like exclusion scheme with nrexcl=1, meaning all atom involved in an
>> bond, angle or dihedral are excluded.
>>
>> You can check the exclusions with
>> $ csg_dump --top XXX --cg YYY --excl
>> and add more bonds to the mapping file to exclude more atom.
>>
>> Christoph
>