create an input file

you will likely want to extract a geometry from a log file first

extract_structure.py filename (on your log file with without the .log extension)

this will give you a file filename.xyz

copy this to filename.inp and add the required infomation (detailed below) to the top of the file

otherwise create a filename.inp file

below is the file for water (water.inp)

info is given on each line

number of atoms

a b c

where

a atom to center the coordinates at (ie will be 0.0 0.0 0.0) here atom 1

b atom you want along the x-axis, in my case atom 2

c atom that will form the xy-plane along with a and b

0.0 0.0 0.0

replace ONE of these with the angle you want to rotate the a-b bond in the xy, xz or yz plane

here I have chosen a 10 degree rotation in the xy plane

currently you can only give ONE non-zero angle here

from then on follows your atomic coordinates (in Angstrom)

note that you can insert a file using vi with the (comand mode):r filename command

 3
 1 2 3
 10.0 0.0 0.0
O   1.995482   0.000000  -0.000000 
H   2.637933  -0.000000   0.697270
H   2.637933   0.000000  -0.697270

run the code

run the code

pos_zmat.py water

- the geom.inp file is an argument (without .inp)

you should see something like the following

[tricia@]/Users/tricia/files/work/jobs/musa $ pos_zmat.py water
directory is /Users/tricia/Files/work/jobs/musa
input file is water.inp
symbolic z-matrix output file is water_osymb.com
xyz coord output file is water_new.xyz
total number of atoms 3
center at atom 1
x-axis along bond to atom 2
atoms in the xy-plane 3
rotation of a-b around xy-pane (z-axis) by 10.0
centered
x-axis aligned
atoms aligned in xy plane
a-b vector rotated
done!

output files

you should see 2 files geom_new.xyz and geom_osymb.xyz

- you can see that atom 1 is at the origin

- you can see that the bond O-H is almost aligned with the x-axis, it has been rotated in the xy plane 10 degrees.

- the molecule is in the xy plane formed by atoms 1,2,3

- note that 3 dummy atoms are added on the position of the x,y and z axies to help with visualisation

geom_new.xyz gives the new coordinate file in the standard format

6   
title 
O    0.000000    0.000000    0.000000 
H    0.933714    0.164639    0.000000 
H   -0.240373    0.917142   -0.000000 
Xx   4.000000    0.000000    0.000000 
Xy   0.000000    4.000000    0.000000 
Xz   0.000000    0.000000    4.000000

geom_osymb.xyz gives the coordinates in the symbolic z-matrix format

%nprocshared=16
%mem=30GB
%chk=water.chk
# opt b3lyp/6-311+g(d,p) int=ultrafine scf=conver=10 empiricaldispersion=gd3bj 
  nosymm geom=z-matrix 
 
 Title Card Required

0 1 
O  0  x1    y1     z1    
H  0  x2    y2     z2    
H  0  x3    y3     z3    
Xx  0  x4    y4     z4    
Xy  0  x5    y5     z5    
Xz  0  x6    y6     z6    

  x1  = 0.000000 
  y1  = 0.000000 
  z1  = 0.000000 
  x2  = 0.933714 
  y2  = 0.164639 
  z2  = 0.000000 
  x3  =-0.240373 
  y3  = 0.917142 
  z3  =-0.000000 
  x4  = 4.000000 
  y4  = 0.000000 
  z4  = 0.000000 
  x5  = 0.000000 
  y5  = 4.000000 
  z5  = 0.000000 
  x6  = 0.000000 
  y6  = 0.000000 
  z6  = 4.000000

pos_zmat.py

#!/opt/local/bin/python3
#
# pos_zmat.py
#
# python script to take cartesian coords and orient a molecule and convert to a symbolic z-matrix format.
# the *.xyz file can be draged into jmol
# the symbolic matrix file is saved in *.com format
#
# to run the script you need an input file *.inp
# if you have a log file run extract_structure.py filename (without log extension) first to get the xyz coordinates
# then copy this file to filename.inp and add the following at the top
# no of atoms
# a b c 
# where
#     a atom to center at
#     b x-axis will align along bond a-b
#     c the atoms a-b-c will form the xy plane
#       this means the  z-axis is perpendicular to the plane of these 3 atoms
# xx yy zz
#     the angle to rotate the a-b bond (vector) in the xy, xz, or yz planes
#     currently only active for one axis at a time!
# then follows a section for atomic symbol/atomic no and coordinates
#
# to run the script type python pos_zmat.py input_file_name (without extension)
#
# this will produce the files
# filename_new.xyz (for jmol)
# filename_osymb.com (for G16)
#
import sys
import os
import math
import numpy as np
dir=os.getcwd()
from numpy import (array, dot, arccos, clip)
from numpy.linalg import norm
#
if len(sys.argv) >= 3:
  print('to run the script please type: python symb_zmat.py input_file_name (without extension)')
  sys.exit()
else:
  base=str(sys.argv[1])
  file=base+'.inp'
  out1=base+'_osymb.com'
  out2=base+'_new.xyz'
  s='directory is '
  print('{0:}{1:}'.format(s,dir))
  s='input file is '
  print('{0:}{1:}'.format(s,file))
  s='symbolic z-matrix output file is '
  print('{0:}{1:}'.format(s,out1))
  s='xyz coord output file is '
  print('{0:}{1:}'.format(s,out2))
#close if
#
# set some defaults
atomtotal=0
atomcenter=0
atomxaxis=0
atomyaxis=0
xa=0.0
xb=0.0
xc=0.0
rotanglex=0.0
rotangley=0.0
rotanglez=0.0
dipole=0.0
nz1=0
nz2=0
nz3=0
a=[]
col_atom=[]
col_x=[]
col_y=[]
col_z=[]
newcol_x=[]
newcol_y=[]
newcol_z=[]
v1=np.zeros(3,dtype=float)
v2=np.zeros(3,dtype=float)
v3=np.zeros(3,dtype=float)

# open file and read lines one by one
f=open(file,"r")

#total number of atoms
line =f.readline()
#print('{0:}'.format(line))
atomtotal=int(line.rstrip())
s='total number of atoms '
print('{0:}{1:}'.format(s,atomtotal))

#read in data center atom, xaxis atom, xy plane atom
line =f.readline()
#print('{0:}'.format(line))
a=line.rstrip().split()
atomcenter=int(a[0])
atomxaxis=int(a[1])
atomyaxis=int(a[2])
s='center at atom '
print('{0:}{1:}'.format(s,atomcenter))
s='x-axis along bond to atom '
print('{0:}{1:}'.format(s,atomxaxis))
s='atoms in the xy-plane '
print('{0:}{1:}'.format(s,atomyaxis))
 
#read in the rotation angles (degrees)
line =f.readline()
b1,c1,d1=line.rstrip().split()
rotplanexy=float(b1)
rotplanexz=float(c1)
rotplaneyz=float(d1)
if rotplanexy != 0:
  s='rotation of a-b around xy-pane (z-axis) by '
  print('{0:}{1:}'.format(s,rotplanexy))
elif rotplanexz != 0:
  s='rotation a-b around xz-pane (y-axis) by '
  print('{0:}{1:}'.format(s,rotplanexz))
elif rotplaneyz != 0:
  s='rotation a-b around yz-pane (x-axis) by '
  print('{0:}{1:}'.format(s,rotplaneyz))
else:
  print('no rotations')
#endif

#read in atoms and coordinates
mx=np.zeros(atomtotal,dtype=float)
my=np.zeros(atomtotal,dtype=float)
mz=np.zeros(atomtotal,dtype=float)
mx1=np.zeros(atomtotal,dtype=float)
my1=np.zeros(atomtotal,dtype=float)
mz1=np.zeros(atomtotal,dtype=float)
mxf=np.zeros(atomtotal,dtype=float)
myf=np.zeros(atomtotal,dtype=float)
mzf=np.zeros(atomtotal,dtype=float)
#print(ny)
n=0
while n < atomtotal :
  line =f.readline()
#  print('{0:}'.format(line))
  a1,b1,c1,d1=line.rstrip().split()
  col_atom.append(a1)
  col_x.append(float(b1))
  col_y.append(float(c1))
  col_z.append(float(d1))
  n=n+1
#close while
#turn list into a numpty array
mx=np.array(col_x)
my=np.array(col_y)
mz=np.array(col_z)
#print(ny)

#center at atomcenter
c=atomcenter-1
#print('center at {0:} {1:9.6f},{2:9.6f},{3:9.6f}'.format(atomcenter,col_x[c],col_y[c],col_z[c]))
mx = mx - mx[c]
my = my - my[c]
mz = mz - mz[c]
print('centered')

#align along x-y atom vector
a=atomcenter-1
b=atomxaxis-1
#print('align along atom vector {0:} {1:}'.format(atomcenter,atomxaxis))
#print('old atom coords {0:} {1:9.4f},{2:9.4f},{3:9.4f}'.format(atomxaxis,col_x[b],col_y[b],col_z[b]))
#print('new atom coords {0:} {1:9.4f},{2:9.4f},{3:9.4f}'.format(atomxaxis,mx[b],my[b],mz[b]))
v0=np.zeros(2,dtype=float)
vxy=np.zeros(2,dtype=float)
vxz=np.zeros(2,dtype=float)
print('x-axis aligned')

#first rotate xy plane
#ie rotate around z-axis
#note that atan2 give the correct quadrant
theta=math.atan2(my[b],mx[b])
t2_deg=math.degrees(theta) 
#print('theta ',round(theta,2), 'radians')
#print('theta ',round(t2_deg,2), 'degrees')
#rotate backwards!
t2_rad=-theta
mx1=mx*np.cos(t2_rad)-my*np.sin(t2_rad)
my1=mx*np.sin(t2_rad)+my*np.cos(t2_rad)
mz1=mz
#print('after xy rotation  {0:9.4f}, {1:9.4f}, {2:9.4f}'.format(mx1[b],my1[b],mz1[b]))

#then rotate in xz plane
#ie rotate "down"around y-axis
theta=math.atan2(mz1[b],mx1[b])
t2_deg=math.degrees(theta) 
#print('theta ',round(theta,2), 'radians')
#print('theta ',round(t2_deg,2), 'degrees')
#rotate backwards!
t2_rad=-theta
mx=mx1*np.cos(t2_rad)-mz1*np.sin(t2_rad)
my=my1
mz=mx1*np.sin(t2_rad)+mz1*np.cos(t2_rad)
#print('after xz rotation  {0:9.4f}, {1:9.4f}, {2:9.4f}'.format(mx[b],my[b],mz[b]))

#then rotate in yz plane
#ie rotate "down"around x-axis
#to bring atom c into the xy plane
b=atomyaxis-1
theta=math.atan2(mz[b],my[b])
t2_deg=math.degrees(theta) 
#print('theta ',round(theta,2), 'radians')
#print('theta ',round(t2_deg,2), 'degrees')
#rotate backwards!
t2_rad=-theta
mx1=mx
my1=my*np.cos(t2_rad)-mz*np.sin(t2_rad)
mz1=my*np.sin(t2_rad)+mz*np.cos(t2_rad)
#print('after yz rotation  {0:9.4f}, {1:9.4f}, {2:9.4f}'.format(mx1[b],my1[b],mz1[b]))
print('atoms aligned in xy plane')

#reset so mx,my,mz contain current values too
mx=mx1
my=my1
mz=mz1

#if rotangle is not zero rotate the bond in xy,xz or yz planes 
t=0
if rotplanexy != 0.0 : 
  a=np.deg2rad(rotplanexy)
  mx1=mx*np.cos(a)-my*np.sin(a)
  my1=mx*np.sin(a)+my*np.cos(a)
  mz1=mz
  t=1
# endif

if t==0 and rotplanexz != 0.0 :
  print('rotate in xz plane')
  a=np.deg2rad(rotplanexz)
  mx1=mx*np.cos(a)-mz*np.sin(a)
  my1=my
  mz1=mx*np.sin(a)+mz*np.cos(a)
  t=1
#endif

if t==0 and  rotplaneyz != 0.0 :
  a=np.deg2rad(rotplaneyz)
  mx1=mx
  my1=my*np.cos(a)-mz*np.sin(a)
  mz1=my*np.sin(a)+mz*np.cos(a)
  t=1
#endif
print('a-b vector rotated')

mfx=mx1
mfy=my1
mfz=mz1

#now write to the _new.xyz file
c=open(out2,"w")
c.write('{0:<3} \n'.format(atomtotal+3))
c.write('{0:<5} \n'.format('title'))
n=0
while n < atomtotal:
  c.write('{0:<3} {1:9.6f}   {2:9.6f}   {3:9.6f} \n'.format(col_atom[n],mfx[n],mfy[n],mfz[n]))
  n=n+1
#close while
#postions dummy atoms on axis
p=4.0
c.write('{0:<3} {1:9.6f}   {2:9.6f}   {3:9.6f} \n'.format('Xx',p,0.0,0.0))
c.write('{0:<3} {1:9.6f}   {2:9.6f}   {3:9.6f} \n'.format('Xy',0.0,p,0.0))
c.write('{0:<3} {1:9.6f}   {2:9.6f}   {3:9.6f} \n'.format('Xz',0.0,0.0,p))
#
#now write to the _symb.com file
c=open(out1,"w")
s="%nprocshared=16\n%mem=30GB\n"
c.write('{0:}'.format(s))
s='%chk='+base+'.chk\n'
c.write('{0:}'.format(s))
s="# opt b3lyp/6-311+g(d,p) int=ultrafine scf=conver=10 empiricaldispersion=gd3bj \n "
c.write('{0:}'.format(s))
s=" geom=z-matrix \n \n "
c.write('{0:}'.format(s))
s="Title Card Required\n\n"
c.write('{0:}'.format(s))
charge=0
multiplicity=1
c.write('{0:} {1:} \n'.format(charge,multiplicity))
#
#write out z-matrix
n=0
while n < atomtotal:
  m=n+1
  c.write('{0:}  0  x{1:<4} y{2:<4}  z{3:<4} \n'.format(col_atom[n],m,m,m))
  n=n+1
#close while
m=n+1
c.write('{0:}  0  x{1:<4} y{2:<4}  z{3:<4} \n'.format('Xx',m,m,m))
m=n+2
c.write('{0:}  0  x{1:<4} y{2:<4}  z{3:<4} \n'.format('Xy',m,m,m))
m=n+3
c.write('{0:}  0  x{1:<4} y{2:<4}  z{3:<4} \n'.format('Xz',m,m,m))
c.write('\n'.format())

#write out variables 
n=0
while n < atomtotal:
  m=n+1
  c.write('  x{0:<3}={1:>9.6f} \n  y{2:<3}={3:9.6f} \n  z{4:<3}={5:9.6f} \n'.format(m,mfx[n],m,mfy[n],m,mfz[n]))
  n=n + 1
#close while
m=n+1
c.write('  x{0:<3}={1:>9.6f} \n  y{2:<3}={3:9.6f} \n  z{4:<3}={5:9.6f} \n'.format(m,p,m,0.0,m,0.0))
m=n+2
c.write('  x{0:<3}={1:>9.6f} \n  y{2:<3}={3:9.6f} \n  z{4:<3}={5:9.6f} \n'.format(m,0.0,m,p,m,0.0))
m=n+3
c.write('  x{0:<3}={1:>9.6f} \n  y{2:<3}={3:9.6f} \n  z{4:<3}={5:9.6f} \n'.format(m,0.0,m,0.0,m,p))
c.write('\n\n'.format())

# close all open files
f.close()
c.close()
print('done!')
exit(0)
#

Mod:Hunt Research Group:position symbolic zmat

Contents

create an input file

run the code

output files

pos_zmat.py

Navigation menu

Personal tools

Namespaces

Variants

Views

More

Search

Navigation

Help

Tools