* imidazole minimization and vibrational analysis
* 

fast on

! output identifier                                                  
set id1 imid_1

open unit 90 write card name @id1.geom

! 1: topology
! 2: parameter
! 3: cutim
! 4: cutnb
! 5: ctonnb
! 6: ctofnb
! 7: electrostatic cutoff regime (shift or switch)
! 8: electrostatic pairing method (atom or group)
! 9: VDW pairing method (vatom or vgroup)
! VDW cutoff regime is vwsitch
set a top_all22_model.inp
set b par_all22_prot.inp
!use infinite cutoffs
set 3 999.0 ! cutim
set 4 999.0 ! cutnb
set 5 800.0 ! ctonnb
set 6 900.0 ! ctofnb
set 7 switch
set 8 atom
set 9 vatom

! read topology and parameters
open unit 9 read form name @a
read rtf card unit 9

open unit 9 read form name @b
read para card unit 9

!generate 
read sequence card
* imidazole
*  
1
imia

generate a first none last none setup warn

!fill ic table with bond lengths and angles from parameter file
ic para

!build molecule
ic seed 1 nd1 1 cg 1 cd2
ic build

coor print

!create nonbond list
update inbfrq -1 ihbfrq 0 -
@7 @8 @9 vfswitch cutnb @4 ctofnb @6 ctonnb @5

mini abnr nstep 50 nprint 10
mini nraph nstep 20 nprint 5

!calculate vibrational spectra, then analyze via MOLVIB
!
vibran
diag
end

!atom numbering
!1  cg
!2  hg
!3  cd2
!4  hd2
!5  nd1
!6  hd1
!7  ce1
!8  he1
!9  ne2

MOLVIB NDI1 21 NDI2 21 NDI3 31 SECO
GFX
PRNT    0
DIM    39   39   57
IC
 1  7  9  0  0   !1  ce1-ne2
 1  9  1  0  0   !2  ne2-cg
 1  1  3  0  0   !3  cg-cd2
 1  3  5  0  0   !4  cd2-nd1
 1  5  7  0  0   !5  nd1-ce1
 1  1  2  0  0   !6  cg-hg
 1  3  4  0  0   !7  cd2-hd2
 1  5  6  0  0   !8  nd1-hd1
 1  7  8  0  0   !9  ce1-he1
 2  5  7  9  0   !10 nd1-ce1-ne2: a1, 5-mem ring
 2  7  9  3  0   !11 ce1-ne2-cd2: a2
 2  9  3  1  0   !12 ne2-cd2-cg: a3
 2  3  1  5  0   !13 cd2-cg-nd1: a4
 2  1  5  7  0   !14 cg-nd1-ce1: a5
 2  3  1  2  0   !15 cd2-cg-hg
 2  5  1  2  0   !16 nd1-cg-hg
 2  9  3  4  0   !17 ne2-cd2-hd2
 2  1  3  4  0   !18 cg-cd2-hd2
 2  9  7  8  0   !19 ne2-ce1-he1
 2  5  7  8  0   !20 nd1-ce1-he1
 2  7  5  6  0   !21 ce1-nd1-hd1
 2  1  5  6  0   !22 cg-nd1-hd1
 3  2  3  5  1   !23 hg-cd2-nd1-cg,   wilson wags to treat hydrogen 
 3  4  1  9  3   !24 hd2-cg-ne2-cd2,  out-of-plane motions
 3  6  1  7  5   !25 hd1-cg-ce1-nd1
 3  8  9  5  7   !26 he1-ne2-nd1-ce1
 4  1  5  7  9   !27 cg-nd1-ce1-ne2, t1
 4  3  1  5  7   !28 cd2-cg-nd1-ne2, t5
 4  9  3  1  5   !29 ne2-cd2-cg-nd1, t4
 4  7  9  3  1   !30 ce1-ne2-cd2-cg, t3
 4  5  7  9  3   !31 nd1-ce1-ne2-cd2, t2
UMAT     0    1    0           ! row normalization
  1  1  1.          2  2  1.          3  3  1.          4  4  1.
  5  5  1.          6  6  1.          7  7  1.          8  8  1.
  9  9  1. 
 10 10  1.         10 11 -0.81       10 12  0.31       10 13  0.31
 10 14 -0.81    
 11 11 -1.12       11 12  1.81       11 13 -1.81       11 14  1.12
 12 27  0.31       12 28 -0.81       12 29  1.0        12 30 -0.81
 12 31  0.31
 13 27 -1.81       13 28  1.12       13 30 -1.12       13 31  1.81
 14 15  1.         14 16 -1.
 15 17  1.         15 18 -1.
 16 19  1.         16 20 -1.
 17 21  1.         17 22 -1.
 18 23  1.
 19 24  1.
 20 25  1.
 21 26  1.
-1
PED      0   15
  1  sCE1NE2     2  sNE2CG      3  sCGCD2      4  sCD2ND1
  5  sND1CE1     6  sHG         7  sHD2        8  sHD1
  9  sHE1       10  dRING      11  dRING'     12  tRING
 13  tRING'     14  dHG        15  dHD2       16  dHE1
 17  dHD1       18  wHG        19  wHD2       20  wHD1
 21  wHE1
 -1
END

stop