NICS计算及其图示|Jerkwin

2014-06-12 14:02:48

使用PQS计算NICS(Nuclear Independent Chemical Shift), 示例输入文件(环丁二烯)如下

TEXT=Cyclobutadiene NICS dense grid for visualization
GEOM=PQS SYMM=0.0
c          -0.793380   -0.661025    0.002932
h           1.560349    1.423031   -0.004224
h          -1.560356   -1.423038    0.004597
h           1.542840   -1.441687   -0.030623
c          -0.785231    0.670514    0.014786
h          -1.542825    1.441702    0.030190
c           0.793371    0.661020   -0.002675
c           0.785233   -0.670517   -0.014981
BASIS=6-311g-dp
OPTImize
SCF DFT=B3LYP
FORCe
JUMP
GEOM PRINT=1
FORCE
HESS
!!!!!!!!!!!!!! NICS Calculation  !!!!!!!!!!!!!!!!
!!                                             !!
!!  For maximum resolution in the final plot,  !!
!!  maintain GRID at 100 and, if the size of   !!
!!  the resulting grid needs to be changed     !!
!!  (ex to cover a larger area), modify STEP.  !!
!!  For instance, the input below produces a   !!
!!  grid of 100x100 NICS points that is 10A    !!
!!  along the sides - where the grid plane     !!
!!  is parallel to the plane of atoms 1, 5,    !!
!!  & 7 but displaced above the plane by 1.0A. !!
!!                                             !!
!!  Note: 100 is the max GRID value            !!
!!                                             !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
NMR NICS=1,5,7 GRID=100 STEP=0.10 DISP=1.0

处理输出文件用于绘图的awk脚本如下

# Language: bash
	File="cyclobutadiene_NICS.out"
	
	awk '' BEGIN{ Npnt=0; Bhr=0.52917720859}
	/GRID for NICS over square/ { Nini=NR; Ndim=$NF
		getline; Size=$(NF-1); sub("=", "", Size); Size=Size*Bhr
	}
	
	Nini!=0 && NR>Nini && /-- Average/ { Npnt++; NICS[Npnt]=-$4 }
	
	END{ for(i=1; i<=Npnt-1; i++) {
		x=int((i-0.1)/Ndim)
		y=i-1-x*Ndim
		printf "%f %f %f\n", x*Size, y*Size, NICS[i]
		}
	}
	' $File >${File%.out}.nics

gnuplot绘图脚本如下

# Language: bash
	gpl set view map;
	set title "NICS of Cyclobutadiene";
	set xlabel "\305"; set ylabel "\305";
	set xtics border in scale 0,0 mirror   norotate offset character 0, 0, 0 autojustify;
	set ytics border in scale 0,0 mirror   norotate offset character 0, 0, 0 autojustify;
	set ztics border in scale 0,0 nomirror norotate offset character 0, 0, 0 autojustify;
	set rtics axis   in scale 0,0 nomirror norotate offset character 0, 0, 0 autojustify;
	set xrange [0:5.25] noreverse nowriteback;
	set yrange [0:5.25] noreverse nowriteback;
	set cbrange [ -20.00000 : 20.00000 ] noreverse nowriteback;
	set cbtics ("Aromatic" -20, "Neutral" 0, "Anti-" 20);
	set palette rgbformulae 33, 13, 10;
	set size square;
	plot "cyclobutadiene_NICS.nics" using 1:2:3 with image title""

环丁二烯