VMD-L Mailing List

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
University of Extremadura
Avda. de Elvas, s/n 06006 Badajoz (Spain)
email: melendez_at_unex.es
phone: +34 924 28 96 55
fax: +34 924 28 96 51