Very primitive code to conflate animated modes; representing modes wi…

…thin modes

JuliaPhonons.jl

@@ -2,7 +2,7 @@ module JuliaPhonons
 
 export atomicmass
 export read_POSCAR, read_meshyaml
-export output_animated_xyz, decompose_eigenmode_atomtype, gnuplot_header, decompose_eigenmode_atom_contributions
+export output_animated_xyz, output_conflated_xyz, decompose_eigenmode_atomtype, gnuplot_header, decompose_eigenmode_atom_contributions
 
 import YAML
 
@@ -118,6 +118,38 @@ function output_animated_xyz(POSCAR::POSCARtype, eigenmode,eigenvector,freq,step
     close(anim)
 end
 
+"Conflated overlapping phonons; a work in proress."
+function output_conflated_xyz(POSCAR::POSCARtype, modecount ,eigenvectors,eigenmodes,steps=200)
+    filename= @sprintf("conflated_%03d.xyz",modecount)
+    anim=open(filename,"w")
+
+    for phi=0:2*pi/steps:2*pi-1e-6 #slightly offset from 2pi so we don't repeat 0=2pi frame
+        phi/=eigenmodes[1] # multiply up to match frequency of first mode
+        # output routines to .xyz format
+        @printf(anim,"%d\n\n",POSCAR.natoms*length(POSCAR.supercell)) # header for .xyz multi part files
+
+        for i=1:POSCAR.natoms
+# Nb: norm of eigenvector is fraction of energy of this mode, therefore you need to 
+#   divide the eigenvector by sqrt(amu) - converting from Energy --> displacement
+            projection=POSCAR.positions[i,:]'
+
+            for (eigenmode,eigenvector) in zip(eigenmodes,eigenvectors)
+                projection+= (1/eigenmode) * eigenvector[i,:]'*sin(phi*eigenmode) / sqrt(atomicmass[POSCAR.atomnames[i]]) # Fractional coordinates
+            end
+
+            projection*=POSCAR.lattice # Scale by lattice [3x3] matrix
+
+            for supercellexpansion in POSCAR.supercell # Runs through which unit cells to print
+                supercellprojection=projection+supercellexpansion'
+                @printf(anim,"%s %f %f %f\n",POSCAR.atomnames[i],supercellprojection[1],supercellprojection[2],supercellprojection[3])
+            end
+        end
+    end
+
+    close(anim)
+end
+
+
 """
 decompose_eigenmode_atomtype(POSCAR::POSCARtype,label,realeigenvector,freq)
 

phonopy_projector.jl

@@ -10,9 +10,15 @@ eigenvectors,eigenmodes=read_meshyaml(open("mesh.yaml"),P)
 myf=open("mode_decomposition.dat","w")
 gnuplot_header(P,f=myf)
 
+output_conflated_xyz(P,1,(eigenvectors[4],eigenvectors[12],eigenvectors[18]),(eigenmodes[4],eigenmodes[12],eigenmodes[18]) ) # generates files anim_{count}.xyz
+
+output_conflated_xyz(P,2,(eigenvectors[4],eigenvectors[7],eigenvectors[10],eigenvectors[13],eigenvectors[16],eigenvectors[19]),(eigenmodes[4],eigenmodes[7],eigenmodes[10],eigenmodes[13],eigenmodes[16],eigenmodes[19]) ) # generates files anim_{count}.xyz
+
+output_conflated_xyz(P,3,(eigenvectors[4],eigenvectors[7],eigenvectors[10],eigenvectors[13],eigenvectors[16],eigenvectors[19],eigenvectors[22],eigenvectors[25],eigenvectors[28],eigenvectors[31],eigenvectors[34]),(eigenmodes[4],eigenmodes[7],eigenmodes[10],eigenmodes[13],eigenmodes[16],eigenmodes[19],eigenmodes[22],eigenmodes[25],eigenmodes[28],eigenmodes[31],eigenmodes[34]) ) # generates files anim_{count}.xyz
+
 for (count,(eigenvector,eigenmode)) in enumerate(zip(eigenvectors,eigenmodes))
     output_animated_xyz(P,count,eigenvector,eigenmode) # generates files anim_{count}.xyz
-    
+
     decompose_eigenmode_atomtype(P,count,eigenvector,eigenmode)
     # Output again, to a file for later printing
     decompose_eigenmode_atomtype(P,count,eigenvector,eigenmode,f=myf)