• Outreach

From: jj mm (jjmm1974_at_yahoo.es)
Date: Thu Dec 09 2010 - 03:11:03 CST

Hi all

I am trying to calculate the phonon spectrum of yttria-doped zirconia. This can be calculated from the velocities of the particles via the VACF of the whole system. I have performed MD simulations using LAMMPS, and proceeded as follows:

1) Generated a list of particles velocities in a multicolumns format ("custom" LAMMPS output)
2) Wrote the velocities list into a .xyz format. After this step I have a file in the standard .xyz format (namely "output.xyz") which contains velocities instead of positions.
3) Calculated the spectrum using the following script:

package require specden

set mol [mol new {output.xyz} waitfor all]

set sel [atomselect \$mol {all}]
set nf [molinfo \$mol get numframes]
set na [\$sel num]

set reslist {}
for {set a 0} {\$a < \$na} {incr a} {
set dlist {}
for {set f 0} {\$f < \$nf} {incr f} {
\$sel frame \$f
lappend dlist [lindex [\$sel get {x y z}] \$a]
}
lassign [specden \$dlist 413.414 3000.0 fourier] flist slist
lappend reslist \$slist
}

# write out the result as: fequency, h1, h2, ...,
set fp [open "dos-1500.dat" "w"]
set ns [llength \$flist]
for {set i 0} {\$i < \$ns} {incr i} {
puts -nonewline \$fp "[lindex \$flist \$i] "
set avg 0.0
for {set a 0} {\$a < \$na} {incr a} {
set val [lindex \$reslist \$a \$i]
#           puts -nonewline \$fp "\$val "
set avg [expr {\$avg + \$val}]
}
puts \$fp "[expr \$avg / \$na]"
}
close \$fp

(I got it from the VMD page. Just commented the "puts..." line because I don't want the DOS associate to each particle, but the global one)

4) After this calculation, I got the file "dos-1500.dat" which, I do not misunderstand, contains the phonon density of states as a function of the wavenumber.

I have a few questions for you guys:

1) Is this procedure right? I mean: it is actually the phonon DOS what I get after the script is run?
2) I made the calculation setting k_max = 3000 cm^-1, and I got 90 frequency points. Is there any way to increase the number of wavenumber values without changing k_max? I am asking this because I calculated the integral below the DOS curve and I found it not to be normalized to 1.
3) Is there any way to do the same thing directly from the LAMMPS trajectory file (i.e., a file which contains positions AND velocities of all the particles of the system)?

Thanks a lot!

Juanjo Melendez
Associate Professor
Department of Physics