Mod:Hunt Research Group/VMDmovie
How to... turn a gaussian optimization into a VMD movie
First create a .xyz file with your coordinates. The format has to be as follows:
<Number of Atoms>
<title>
<Symbol or Number of Atom1> <x> <y> <z>
<Symbol or Number of Atom2> <x> <y> <z>
<Symbol or Number of Atom3> <x> <y> <z>
...
repeat this for every frame of your animation (= every step of your optimization). You cannot change the number of atoms or the title after the first step, but you still have to include them (or at least a blank line for each).
To simplify this process you can get a copy of the XYZmovie.py script. This will do a number of things:
1.) It will read every geometry from a gaussian output file and bring it into the correct XYZ format for VMD. If you do this with a frequency or singlepoint calculation you are just going to get a single geometry, if you try it with an optimization, you will get a series of frames, one for each step of the optimization.
2.) It will write the output into a file. :)
3.) It will move one atom of your choice to the origin of the coordinate system.
4.) It will rotate the molecule so that a second atom of your choice is on the y axis.
5.) It will rotate the molecule so that a third atom of your choice is in the y,z plane.
6.) It will rotate the molecule so that the second and third atom have a positive y value.
Point 3 to 6 are in order to stop the molecule from 'hopping around'. The atoms should be chosen so that they are not strongly affected by the movement in your molecule, otherwise the animations may look 'odd'.
The syntax for the script is
XYZmovie.py <input file> <output file> <atom 1> <atom2> <atom3>
The atoms are counted starting from 0!
The resulting file can then be opened in VMD and should result in a trajectory/movie. The representation can then be changed as usual in VMD. To save the result as a movie, open Extensions -> Visualization -> Movie Maker in VMD.
In the new window the following options should be set:
Renderer should be set to either 'Snapshot' or 'Tachyon (Raytracer)'. The later gives slightly better looking pictures but takes significantly longer to render. If time is any issue at all, go for Snapshot.
Movie Settings: check 'Trajectory' and uncheck 'Delete image files'
Format should be 'Animated GIF (ImageMagick)'
Set the working directory and give the whole thing a name of your choice. Finally click on Make Movie.
This will create a range of files of the patter <directory>/<name>.<frame number>.ppm, each a picture of another frame, as well as the file <directory>/<name>.gif with the video. If this video is to fast or slow, you can create a new video from the frames by running
convert -delay
from a console window. The standard delay is 4.17. You can copy the command from the VMD console window to make it easier. A delay of >10 will most likely result in a video which starts to 'jump', not showing a smooth motion.
P.S.: If the optimization itself is jerking (i.e. goes back and forth, 'shivers' etc.) in the VMD window before exporting it as a video, open Graphics -> Representations -> Trajectory and increase the Trajectory Smoothing Window Size. A value of 5 or 10 worked well for me, but try this for yourself. The higher the value, the more detail you are going to 'drown', up to a point where the rotation of a methyl group will actually move the hydrogens through the carbon.
Python script
#!/usr/bin/python
import sys, re, math
def shift(matrix,index=0):
x=matrix[index][2]
y=matrix[index][3]
z=matrix[index][4]
for i in range(len(matrix)):
matrix[i][2] = matrix[i][2] - x
matrix[i][3] = matrix[i][3] - y
matrix[i][4] = matrix[i][4] - z
return matrix
def flat(matrix,index,ix,iy,iz):
matrix=makekugel(matrix,ix,iy,iz)
y=matrix[index][3]
for i in range(len(matrix)):
matrix[i][3] = matrix[i][3] - y
matrix=makekart(matrix,ix,iy,iz)
return matrix
def kxflat(matrix,index=1,axis=3,rota=2,rotb=4):
winkel=-math.atan(matrix[index][rota]/matrix[index][rotb])
for i in range(len(matrix)):
x=matrix[i][2]
y=matrix[i][3]
z=matrix[i][4]
matrix[i][2] = x * math.cos(winkel) + z * math.sin(winkel)
matrix[i][3] = y
matrix[i][4] = z * math.cos(winkel) - x * math.sin(winkel)
if matrix[index][4] < 0 :
for i in range(len(matrix)):
matrix[i][2] = - matrix[i][2]
matrix[i][3] = matrix[i][3]
matrix[i][4] = -matrix[i][4]
return matrix
def kzflat(matrix,index=1,axis=2,rota=4,rotb=3):
winkel=-math.atan(matrix[index][rota]/matrix[index][rotb])
for i in range(len(matrix)):
x=matrix[i][2]
y=matrix[i][3]
z=matrix[i][4]
matrix[i][2] = x
matrix[i][3] = y * math.cos(winkel) - z * math.sin(winkel)
matrix[i][4] = z * math.cos(winkel) + y * math.sin(winkel)
return matrix
def makekugel(matrix,ix,iy,iz):
for i in range(len(matrix)):
x = math.sqrt(matrix[i][ix]**2+matrix[i][iy]**2+matrix[i][iz]**2)
if x==0:
y=0
z=0
else:
y = math.atan2(matrix[i][iy],matrix[i][ix])
z = math.acos(matrix[i][iz]/(math.sqrt(matrix[i][ix]**2+matrix[i][iy]**2+matrix[i][iz]**2)))
matrix[i][2] = x
matrix[i][3] = y
matrix[i][4] = z
return matrix
def makekart(matrix,ix,iy,iz):
for i in range(len(matrix)):
x = matrix[i][2] * math.sin(matrix[i][4]) * math.cos(matrix[i][3])
y = matrix[i][2] * math.sin(matrix[i][4]) * math.sin(matrix[i][3])
z = matrix[i][2] * math.cos(matrix[i][4])
matrix[i][ix] = x
matrix[i][iy] = y
matrix[i][iz] = z
return matrix
def geomread():
global atoms
atoms=[]
#mull=open(file,'r')
#xyz=open(outfile,'w')
while True:
line=infile.readline()
if 'Standard orientation' in line:
#print("hurra")
atoms=[]
for i in range(0,5):
line=infile.readline()
#print(line)
while True:
if '--' in line:
break
items=line.split()
atoms.append([int(items[1]),float(items[1]),float(items[3]),float(items[4]),float(items[5])])
line=infile.readline()
#print(atoms)
atoms=shift(atoms,int(sys.argv[3]))
atoms=kxflat(atoms,int(sys.argv[4]))
atoms=kzflat(atoms,int(sys.argv[4]))
atoms=kxflat(atoms,int(sys.argv[5]))
outfile.write(str(len(atoms))+"\n")
outfile.write(sys.argv[1]+"\n")
for i in atoms:
outfile.write(str(i[0])+" "+str(i[2])+" "+str(i[3])+" "+str(i[4])+"\n")
if line=='':
break
return()
infile = open(sys.argv[1],'r')
outfile = open(sys.argv[2],'w')
print sys.argv[1]
print sys.argv[2]
#print(atoms)
#print(charge)
if sys.argv[1].endswith(".log"):
geomread()
#quit()