0

I have a Lammps Dumpfile in the format "id mol type xu yu zu", where id is the ID of each molecule and mol is the ID of each molecule (atoms connected by a bond).

import MDAnalysis as mda
from MDAnalysis.transformations import boxdimensions

dumpfile_name = "dumpfile.dat"
dim = [20, 20, 20, 90, 90, 90]#Lx,Ly,Lz,alpha,beta,gamma

u = mda.Universe(dumpfile_name, atom_style='id mol type xu yu zu',topology_format='LAMMPSDUMP')
transform = boxdimensions.set_dimensions(dim)
u.trajectory.add_transformations(transform)

At this point I'd like to do an analysis for each molecule, e.g. the radius of gyration and let's say I have 4 of these molecules and I want to calculate the average over time. What is the correct syntax to achieve something like in the following (5th line).

N_Molecules = 4
for ts in u.trajectory:
   Rgyr = 0
   for i in range(1,N_Molecules+1): #LAMMPS counts IDS starting from 1
      molecule = u.select_atoms(f"mol {i}") #same syntax as chosing e.g. type
      Rgyr += molecule.radius_of_gyration()
   Rgyr /= N_Molecules
   print(f"Rgyr {Rgyr:.2f}")

Here a small test dumpfile "dumpfile.dat":

ITEM: TIMESTEP
200000
ITEM: NUMBER OF ATOMS
20
ITEM: BOX BOUNDS pp ff ff
0.0000000000000000e+00 2.0000000000000000e+01
0.0000000000000000e+00 2.0000000000000000e+01
0.0000000000000000e+00 2.0000000000000000e+01
ITEM: ATOMS id mol type xu yu zu
10 2 1 1.80692 16.9556 2.98775
6 2 1 3.54256 17.9352 3.37097
9 2 1 1.27783 17.3512 3.67999
8 2 1 2.05019 17.7925 4.11131
7 2 1 2.96428 17.6825 4.12718
1 1 1 1.77296 1.8086 4.67644
3 1 1 1.59248 3.26525 4.73281
2 1 1 1.74302 2.48995 5.3763
4 1 1 1.61248 3.73294 5.57016
5 1 1 2.50437 4.02051 5.94323
12 3 1 10.2844 1.90882 8.50567
11 3 1 10.1476 2.67261 7.88784
13 3 1 10.6363 1.62 9.32461
15 3 1 11.7797 0.655019 10.2963
14 3 1 11.1027 1.33187 10.1537
18 4 1 7.72979 17.8708 15.0419
17 4 1 7.87695 18.7577 15.3529
20 4 1 7.24772 17.0584 16.0629
19 4 1 8.03517 17.4917 15.8818
16 4 1 7.17389 19.1475 15.8897
ITEM: TIMESTEP
300000
ITEM: NUMBER OF ATOMS
20
ITEM: BOX BOUNDS pp ff ff
0.0000000000000000e+00 2.0000000000000000e+01
0.0000000000000000e+00 2.0000000000000000e+01
0.0000000000000000e+00 2.0000000000000000e+01
ITEM: ATOMS id mol type xu yu zu
13 3 1 1.46292 6.29151 4.25376
15 3 1 0.280798 7.44181 3.77524
12 3 1 1.17366 5.4189 4.68029
11 3 1 2.1382 5.50442 4.85405
14 3 1 0.605659 6.78431 4.37785
7 2 1 1.65536 17.7064 8.32545
1 1 1 2.15224 3.62475 8.7687
10 2 1 2.44988 16.1844 9.11292
8 2 1 1.46369 16.8784 8.83759
6 2 1 2.43784 17.764 8.821
3 1 1 1.91573 3.01187 10.0432
2 1 1 2.41498 3.78368 9.68284
9 2 1 2.02099 16.859 9.59798
5 1 1 0.274499 2.43497 10.4381
4 1 1 1.22946 2.64076 10.5169
17 4 1 8.02326 17.7633 13.6459
16 4 1 7.38805 18.4861 13.6313
18 4 1 8.4123 18.6222 13.3534
20 4 1 7.64567 19.4217 13.3507
19 4 1 8.09438 19.0231 14.1376

This leads to the following error:

MDAnalysis.exceptions.SelectionError: Unknown selection token: 'mol'

Using the token "id"/ "type" works. Now I could write my own method to get the mol_id from the id token, but this feels like not being the intended way.

ludwig_boltzplatz
  • 21
  • 3
  • I notice that you define a variable named `dumpfilename` but have `dumpfile_name` as an argument to `mda.Universe()`. Typo? – jjramsey Aug 24 '22 at 16:41
  • Also, what is `mon` supposed to be? – jjramsey Aug 24 '22 at 16:42
  • You're right I tried to write a condensed example and introduced these errors. I edited the question. – ludwig_boltzplatz Aug 24 '22 at 16:47
  • Even with the typos fixed, it's still hard to tell what you are asking for. Offhand, there isn't any *syntax* problem that jumps out at me, and "something like in the following" is just vague. How is your code not already doing what you want? – jjramsey Aug 24 '22 at 17:15
  • I made some changes to the original post so that it is a working example now. – ludwig_boltzplatz Aug 25 '22 at 09:37
  • So the syntax you're asking about is the syntax of the argument for the `select_atoms()` method? – jjramsey Aug 25 '22 at 13:16

0 Answers0