[CP2K-user] How to plot RDFs for NPT simulations?

[Travis]

Hi,

There is nothing unique about an rdf from NPT. The box size is needed only 
to unwrap everything correctly. Use DCD for MD *exclusively *in the future. 
PDB is also an option but it's rather low precision. You can convert using 
Atomic Simulation Environment to a format readable by VMD that allows rdf 
with variable cells.

My directory contains the files, 
foo.cell  foo-pos-1.xyz

I generated random unit cell parameters with,
for ((i=0; i<310; i++)); do a=`echo "30 + $RANDOM/30" | bc`; b=`echo "30 + 
$RANDOM/30" | bc`; c=`echo "30 + $RANDOM/30" | bc`; echo "$i $i $a $b $c 90 
90 90"; done > foo.cell

In bash, $RANDOM is any 16-bit integer (ceil is 32767); I add 30 just in 
case $RANDOM is zero. I don't remember the exact format of the .cell file, 
so I took a guess of it as `step time, a, b, c, alpha, beta, gamma`.

Here are some of the lines,
0 0 366 530 687 90 90 90
1 1 146 367 536 90 90 90
2 2 928 744 656 90 90 90
3 3 767 1043 943 90 90 90
4 4 621 932 711 90 90 90
5 5 854 1113 581 90 90 90
6 6 1033 204 505 90 90 90
7 7 1057 179 1028 90 90 90
8 8 222 944 814 90 90 90
9 9 552 339 296 90 90 90
10 10 381 807 488 90 90 90
11 11 1020 745 763 90 90 90
12 12 1121 126 34 90 90 90
13 13 802 552 413 90 90 90
14 14 388 1115 473 90 90 90
15 15 752 171 48 90 90 90
16 16 963 438 595 90 90 90
17 17 1039 374 766 90 90 90
18 18 858 514 198 90 90 90
19 19 970 710 161 90 90 90
20 20 175 505 992 90 90 90

Format doesn't matter, just needs to be in columns (with 2 buffer columns 
with step and physical time listed).

And now a Python script for use with ASE,
#!/usr/bin/env python3

'''
usage: python3 xyzcell2gro.py

Reads XMOL format file, adds lattice parameters, dumps Gromos file that can
be read by VMD.
'''

import glob
import os
import sys
import numpy as np

import ase
from ase.io import read, write

dir = os.getcwd()

# trajectory is XMOL format
# cell       is step time a b c alpha beta gamma
for file in glob.glob('%s/*-pos-*.xyz' % dir):
   filename, fileext = os.path.splitext(file)
   fileroot, filescrap = filename.split('-pos')
   trj = read('%s%s' % (filename, fileext) , format='xyz', index=':')
   cel = np.loadtxt('%s.cell' % (fileroot), usecols=range(2,8), dtype=np.
double)
   nframes = len(trj)
   for frame in range(nframes):
      trj[frame].set_cell([ cel[frame][0], cel[frame][1], cel[frame][2], cel
[frame][3], cel[frame][4], cel[frame][5] ])
      trj[frame].set_pbc([True, True, True])
      write('%s.g96' % (fileroot), trj[frame], format='gromos', append=True)

The gromos format needs extension .g96 to be read by default in VMD.

-T


On Wednesday, October 23, 2019 at 2:27:46 PM UTC-3, Ant wrote:
>
> If anyone knows how to make RDFs for NPT simulation results that already 
> exist, it would be great to know.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20191023/3ffd6dd0/attachment.htm>