Talk:Mod:Hunt Research Group/lam fold xyz

From ChemWiki
Jump to: navigation, search

A script to convert a LAMMPS traj file to an xyz file and fold and center the cell

  • copy and paste the script into a file called lam_fold.py
  • run by typing: python lam_fold.py lammps_file XYZ_file no_steps
lammps_file is the lammps trajectory file
XYZ_file is the output xyz file
no_steps is the number of steps to cut from a large file
this will also produce a file data_out that contains the atom names used
==A script to convert a LAMMPS traj file to an xyz file==
*copy and paste the script into a file called lam2xyz.py

*you will need to alter the atomList to match the 2letter names for your atoms!

*run by typing: '''python lam_cut.py lammps_file XYZ_file no_steps'''
::lammps_file is the lammps trajectory file
::XYZ_file is the output xyz file
::no_steps is the number of steps to cut from a large file

<pre>
#!/usr/bin/env python
import sys,string

# run this script by typing python lam_fold lammps_traj_file *.xyz no_steps 
# lammps_file is the input lammps trajectory file
# *.xyz is the output xyz coordinates file
# no_steps is the number of trajectory frames to record
# will also produce a data file with the atom label list

# arguments from running the file,
# argv[1]=lammps file
# argv[2]=output *.xyz file
# argv[3]=number of frames to cut
inputFile=open(sys.argv[1],'r')
outFile=open(sys.argv[2],'w')
no_frames=int(sys.argv[3])
dataFile=open('data_out','w')

#lammps file foromat is
#title time step followed by timestep
#  exammple:
#  ITEM: TIMESTEP
#  0
#title number of atoms folled by the number
# exammple:
# ITEM: NUMBER OF ATOMS
# 14608
#title box information followed by numbers
#where xlo and xhi define the box dimension ends (in this case for x axis)
#  example:
#  ITEM: BOX BOUNDS pp pp pp
#  1.8280820780474468e+01 6.1719179219536784e+01
#  2.7421231170696949e+01 9.2578768829227556e+01
#  2.5136128573138141e+01 8.4863871426850039e+01
#then the list of atoms and their coordinates
#  example:
#  ITEM: ATOMS element xu yu zu
#  N 24.4607 26.7811 30.8739
#  O 25.8947 28.7289 31.7749
#  O 22.0142 26.5087 30.1196

#the VMD xyz file format is
#number of atoms
#comment (normally the time step)
#atom_name x y z

#now lets start the conversion
#loop through reading the lines
#if this line has content do the following
#otherwise quit
count=int(1)
atomLabels = [ ]
print('atomLabels {0}'.format(atomLabels))
print('no_frames {0}'.format(no_frames))
while count <= no_frames:
   print('count {0}'.format(count))

#  timestep
   line=inputFile.readline()
   line=inputFile.readline()
   line=line.rstrip().split()
   time_step=line[0]
   print('timestep is {0}'.format(time_step))

#  number of atoms
   line=inputFile.readline()
   line=inputFile.readline()
   line=line.rstrip().split()
   total_atoms=int(line[0])
#   print('total atoms is {0}'.format(total_atoms))

#  then write data to the *.xyz file
   outFile.write('{0}{1}'.format(total_atoms,'\n'))
   outFile.write('time step {0}{1}'.format(time_step,'\n'))

#  box dimensions
   line=inputFile.readline()
   line=inputFile.readline()
   line=line.rstrip().split()
   x_min=float(line[0])
   x_max=float(line[1])
   line=inputFile.readline()
   line=line.rstrip().split()
   y_min=float(line[0])
   y_max=float(line[1])
   line=inputFile.readline()
   line=line.rstrip().split()
   z_min=float(line[0])
   z_max=float(line[1])
   print('x_min={0:.2f}, x_max={1:.2f}'.format(x_min,x_max))
   print('y_min={0:.2f}, y_max={1:.2f}'.format(y_min,y_max))
   print('z_min={0:.2f}, z_max={1:.2f}'.format(z_min,z_max))
   x_side=x_max-x_min
   y_side=y_max-y_min
   z_side=z_max-z_min
   print('x_side={0:.2f}, y_side={1:.2f}, z_side={2:.2f}'.format(x_side,y_side,z_side))
   x_center=x_min + x_side/2.0
   y_center=y_min + y_side/2.0
   z_center=z_min + z_side/2.0
   print('x_center={0:.2f}, y_center={1:.2f}, z_center={2:.2f}'.format(x_center,y_center,z_center))

#   reading atom name and coordinates
#   write out atom name and xyz coorinates 
   line=inputFile.readline()
   i=1
   while i <= total_atoms:
      line=inputFile.readline()
      line=line.rstrip().split()
      atom_name=line[0]
      x_atom=float(line[1])
      y_atom=float(line[2])
      z_atom=float(line[3])
#     adjusting if outside cell dimensions
      if x_atom > x_max : x_atom = x_atom - x_max
      if x_atom < x_min : x_atom = x_max - x_atom
      if y_atom > y_max : y_atom = y_atom - y_max
      if y_atom < y_min : y_atom = y_max - y_atom
      if z_atom > z_max : z_atom = z_atom - z_max
      if z_atom < z_min : z_atom = z_max - z_atom
#     adjusting cell center to origin
      x_atom = x_atom - x_center
      y_atom = y_atom - y_center
      z_atom = z_atom - z_center
      outFile.write('{0:2s}  {1:8.4f}  {2:8.4f}  {3:8.4f}  {4}'.format(atom_name,x_atom,y_atom,z_atom,'\n'))
#
      if count == 1:
         if not line[0] in atomLabels :
           atomLabels.append(line[0])
#           print('added {0} to atomLabels'.format(line[0]))
           dataFile.write('{0}  '.format(line[0]))
      i=i + 1
   count=count + 1
#endwhile

#close the open files
inputFile.close()
outFile.close()
dataFile.close()
#print('requested number of frames = {0}'.format(no_frames))
print('printed number of frames = {0}'.format(count-1))
print('atom labels can be found in data_out {0}'.format(atomLabels))