Mod:Hunt Research Group:molecules in an electric field

Molecules in an electric field

introduction

getting started with the field keyword

which optimisation to use

accelerating optimisation

using the field keyword

• we apply a field along an axial direction, x or y or z

field=x+250

-for field=x+10 a field of 0.001au is applied

-for more options goto the manual page here link

• units

-conversion for an electric field: 1 au = 5.14220652 × 10^11 V/m =51.4 V/Å =514 V/nm

-so 0.001au = 51.4/1000 V/Å = 0.0514 V/Å or 0.514V/nm

-this is equivalent to 0.05V/Å

-to apply 1V/nm we need approximately a field=x+20 (2*0.514V/nm=1.028V/nm)

optimising in a field

• the field is applied along an axis direction

-use the show axis option in gaussview to show this direction

-by convention electric field arrows point from postive to negative

-a positive x-axis field=+250 for example field will repel a negative ion (move it to -x)

• the molecule must not reorient when setting up the job

-position the molecule explicitly where you want it within the coordinated system

-you MUST use nosymm

-things will go faster if you orient the molecular dipole with the field

-in chemistry the dipole points positive to negative (in physics it can be the reverse!)

-to get the correct dipole vector you must optimise first under the same conditions as will be applied with the field

• we have to be very careful in which cooridnates we optimise

-this is because internal coordinates remove the center of mass motion which is relevant in an electric field

-you MUST use opt=z-matrix

-it is better to use symbolic cartesian coordinates

-below is an example for acetate

-the first zero tells gaussian that symbolic cartesian coordinates follow

-if you run a charged molecule the molecule will move along the electric field

-to avoid this fix one atom in the molecule to the center of the axial system

-in the z-matrix format variables (to be optimised) are listed first and fixed constant are listed second

-so the coordinates of the last atom below are held constant (at the origin)

%chk=ch3co2.chk
# opt=(z-matrix,maxcycles=100) hf/3-21g nosymm field=x+250

Title Card Required

-1 1
6   0   x1 y1 z1
8   0   x2 y2 z2 
8   0   x3 y3 z3 
6   0   x4 y4 z4 
1   0   x5 y5 z5 
1   0   x6 y6 z6 
1   0   x7 y7 z7 

x2=-0.544205
y2=-1.144935
z2= 0.029999
x3=-0.568491
y3= 1.132502 
z3= 0.024945
x4= 1.542020
y4= 0.011495
z4=-0.117920
x5= 1.805270
y5=-0.088033
z5=-1.168772
x6= 1.945393
y6= 0.945207
z6= 0.252204
x7= 1.965502
y7=-0.830194
z7= 0.416379

x1= 0.000000
y1= 0.000000 
z1= 0.000000

• refence state

-make sure you run a reference state calculation without the field

-this has the bonus of giving the dipole moment as well

-currently we cannot orientate the dipole, but hopefully soon

-first optimise the molecule (with the z-matrix as defined for the field calculation)

# opt=(z-matrix,maxcycles=100) hf/3-21g nosymm

a code to orient the molecule

• copy the code pos_zmat into a file called pos_zmat.py in your bin directory

pos_zmat.py link

• change the first line shebang to direct to YOUR python3

mine is on my mac #!/opt/local/bin/python3

on raapoi is #!/usr/bin/env python3.8

• extract the xyz coordinates for your molecule

use "extract_structure.py"

that is the xyz coordinates in the form atom_symbol x y z

for example here are the coordinates of a water molecule

note that we require atomic symbols not numbers

O   1.995482   0.000000  -0.000000 
H   2.637933  -0.000000   0.697270 
H   2.637933   0.000000  -0.697270

• create a file called water.inp

below is the file for water

info is given on each line

1. number of atoms

2. atom_1 atom_2 atom_3

3. angle_1 angle_2 angle_3

atom_1 =atom to center the coordinates at (ie will be 0.0 0.0 0.0) in my case I want to center at atom 1

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

atom_3 =atom to form the xy plane through atom_1, atom_2 and atom_3, the z-axis will be perpendicular to this plane

currently you can only give ONE angle the others must be set to 0.0

angle_1 rotatare in xy

angle_2 rotate in xz

angle_3 rotate in yx

4. and on are your coordinates

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

the water.inp file as an argument (without .inp)

pos_zmat.py water

you should see something like the following

directory is /Volumes/Tricia_Home/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
write .xyz file
done!

• output files

you should see 2 files water_new.xyz and water_osymb.xyz

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

O    0.000000    0.000000    0.000000 
H    0.915812    0.245391    0.000000 
H   -0.074822   -0.020048   -0.944949 
Xx   4.000000    0.000000    0.000000 
Xy   0.000000    4.000000    0.000000 
Xz   0.000000    0.000000    4.000000 

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

you will add these into your com file

you can see that atom 1 is at the origin

you can see that atom 2 has coordinates in the x and y direction very close to zero

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

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.915812 
  y2  = 0.245391 
  z2  = 0.000000 
  x3  =-0.074822 
  y3  =-0.020048 
  z3  =-0.944949 
  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 

• create a com file

use the vi insert file option (comand mode):r geom_osymb.xyz

here is my new water_field.com file

you MUST visualise this in gaussview to make sure everything is ok!

%nprocshared=3
%mem=3GB
%chk=water_osymb.chk
#b3lyp/3-21G int=ultrafine scf=conver=10
 opt=(maxcyc=50)  nosymm field=x+250

 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
:

and here is my gaussview check

with the x-axis to atom 2 angle of 10 degrees

scan a coordinate

• you need to build a z-matrix file that includes symbolic coordinates and internal coordinates together
• for z-matrix the first set of parameters are optimised and the second set are frozen
• shown by the example for water below

atom1 is fixed in place with symbolic coords (note the 0 after the O label)

O-H distance is allowed to vary (r1)

the angle H-O-H is opened from 100º in 5 steps by 15º each time

the scanned coordinate is "frozen" on each step so goes into the block of frozen variables

%nprocshared=5
%mem=3GB
%chk=water_scan_osymb.chk
#b3lyp/3-21G nosymm field=x+50
 opt=(z-matrix,maxcyc=5)

 Title Card Required

0 1
O  0  x1    y1     z1
H  1  r1
H  1  r1  2 a1

  r1  1.0
  
  x1  = 0.000000
  y1  = 0.000000
  z1  = 0.000000
  a1  100.0 S 2 15.0

• a torsion angle rotation is more difficult
• unfortunately you cannot read the z-matrix from a checkpoint file so you need to start from a pre-optimised geometry and then set up the z-matrix rotation
• here is a demo job for alkyl rotation of emim

the C2 is set at the origin (0,0,0) and the coords are frozen to keep the coordinate center fixed

the ring atoms are also all frozen in cartesian coords

then we need to define all the "moving" atoms with internal coordinates

lastly I have the dummy atoms which are showing were the axes are

notice the separation of active and frozen variables

the scan coordinate is given LAST we are rotating 18 steps of 10 degrees

• obviously for your application you will want to improve the method/basis-set and maxcycles!

%nprocshared=5
%mem=5GB
%chk=emim_outplane_scan2_osymb.chk
#b3lyp/3-21G nosymm field=x+50
 opt=(z-matrix,maxcycles=25)

 Title Card Required

1 1
6  0  x1    y1     z1
1  0  x2    y2     z2
6  0  x3    y3     z3
1  0  x4    y4     z4
6  0  x5    y5     z5
1  0  x6    y6     z6
7  0  x7    y7     z7
6  0  x8    y8     z8
1  0  x9    y9     z9
1  0  x10   y10    z10
1  0  x11   y11    z11
7  0  x12   y12    z12
6  12 rnc   5   acnc   3  dcncc
1  13 rch1  12  ahcn1  1  dhcn1
1  13 rch2  12  ahcn2  1  dhcn2
6  13 rcc   12  accn   1  dccnc
1  16 rch3  13  ahcc1  12 dhccn1
1  16 rch4  13  ahcc2  12 dhccn2
1  16 rch5  13  ahcc3  12 dhccn3
Xx  0  x20   y20    z20
Xy  0  x21   y21    z21
Xz  0  x22   y22    z22

  rnc    1.481
  acnc   126.0
  dcncc  180.0
  rch1   1.09
  rch2   1.09
  ahcn1  107.0
  ahcn2  107.0
  dhcn1  -135.0
  dhcn2  -20.0
  rcc    1.524
  accn   112.0
  rch3   1.09
  rch4   1.09
  rch5   1.09
  ahcc1 112.0
  ahcc2 112.0
  ahcc3 112.0
  dhccn1 -60.0
  dhccn2 180.0
  dhccn3 60.0

  x2  = 1.076881
  y2  = 0.000000
  z2  = 0.000000
  x3  =-2.100705
  y3  = 0.680775
  z3  = 0.000000
  x4  =-2.914688
  y4  = 1.384929
  z4  =-0.004197
  x5  =-2.099251
  y5  =-0.680237
  z5  =-0.004583
  x6  =-2.913526
  y6  =-1.383789
  z6  =-0.010663
  x7  =-0.779163
  y7  = 1.086932
  z7  = 0.003562
  x8  =-0.313764
  y8  = 2.480517
  z8  =-0.001050
  x9  = 0.773736
  y9  = 2.489310
  z9  = 0.017415
  x10 =-0.666276
  y10 = 2.977026
  z10 =-0.904172
  x11 =-0.696108
  y11 = 2.991482
  z11 = 0.881532
  x12 =-0.778234
  y12 =-1.085961
  z12 =-0.008125
  x1  = 0.000000
  y1  = 0.000000
  z1  = 0.000000
  x20 = 4.000000
  y20 = 0.000000
  z20 = 0.000000
  x21 = 0.000000
  y21 = 4.000000
  z21 = 0.000000
  x22 = 0.000000
  y22 = 0.000000
  z22 = 4.000000
  dccnc 100.0 S 18 10.0