* imidazole minimization and geometric 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

! autogenerate the ic table and put into topology (reformat as required)
!ic generate

!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

!assign atom numbers to variable for quick commands
coor stat sele type cg end
set cg ?selatom
coor stat sele type hg end
set hg ?selatom
coor stat sele type cd2 end
set cd2 ?selatom
coor stat sele type hd2 end
set hd2 ?selatom
coor stat sele type nd1 end
set nd1 ?selatom
coor stat sele type hd1 end
set hd1 ?selatom
coor stat sele type ce1 end
set ce1 ?selatom
coor stat sele type he1 end
set he1 ?selatom
coor stat sele type ne2 end
set ne2 ?selatom

!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

coor print

!variables for statistics on bonds
set bdif 0.0
set bdif2 0.0
set count 0

write title unit 90
* Bond section
* atoms, charmm, target, difference
*

!note that the experimental distances have been
!corrected for the experimental data being based
!on the NE2 atom being protonated versus ND1 for
!the CHARMM convention.

!experimental data based on coordinates in Table 5
!of Christen et al., Z.Naturforsch 37a:1378-1385 (1982)
!ce1-ne2

!statistical analysis includes hydrogens; often the
!nonhydrogen analysis is done separately from the hydrogens.

set exp 1.314
quick @ce1 @ne2
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ce1-ne2: ?dist @exp @dif
*

!ce1-nd1
set exp 1.364
quick @ce1 @nd1
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ce1-nd1: ?dist @exp @dif
*

!ne2-cd2
set exp 1.382
quick @cd2 @ne2
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ne2-cd2: ?dist @exp @dif
*

!nd1-cg
set exp 1.377
quick @nd1 @cg
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ne2-cd2: ?dist @exp @dif
*

!cg-cd2
set exp 1.364
quick @cg @cd2
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*cg-cd2: ?dist @exp @dif
*

!ce1-he1
set exp 1.079
quick @ce1 @he1
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ce1-he1: ?dist @exp @dif
*

!nd1-hd1
set exp 0.998
quick @nd1 @hd1
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*nd1-hd1: ?dist @exp @dif
*

!cg-hg
set exp 1.078
quick @cg @hg
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*cg-hg: ?dist @exp @dif
*

!cd2-hd2
set exp 1.077
quick @cd2 @hd2
calc dif = ?dist - @exp
calc bdif = @bdif + abs ( @dif )
calc bdif2 = @bdif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*cd2-hd2: ?dist @exp @dif
*

!calculate statistics
! average difference
calc avediff = @bdif / @count
calc avediff2 = @bdif2 / @count
calc stdev = sqrt ( ( @avediff2 - ( @avediff * @avediff ) ) / (@count - 1) )

write title unit 90
* bond difference statistics
* average difference: @avediff
* standard deviation: @stdev
* ____________________________
*

!variables for statistics on bonds
set adif 0.0
set adif2 0.0
set count 0

write title unit 90
* Angle section
* atoms,    charmm, target, difference
*

!ce1-ne2-cd2
set exp 104.94
quick @ce1 @ne2 @cd2
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ce1-ne2-cd2: ?thet @exp @dif
*

!ne2-cd2-cg
set exp 110.69
quick @ne2 @cd2 @cg
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ne2-cd2-cg: ?thet @exp @dif
*

!cd2-cg-nd1
set exp 105.48
quick @cd2 @cg @nd1
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*cd2-cg-nd1: ?thet @exp @dif
*


!cg-nd1-ce1
set exp 106.90
quick @cg @nd1 @ce1
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*cg-nd1-ce1: ?thet @exp @dif
*

!nd1-ce1-ne2
set exp 111.99
quick @nd1 @ce1 @ne2
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*nd1-ce1-ne2: ?thet @exp @dif
*

!nd1-ce1-he1
set exp 122.51
quick @nd1 @ce1 @he1
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*nd1-ce1-he1: ?thet @exp @dif
*

!ne2-ce1-he1
set exp 125.50
quick @ne2 @ce1 @he1
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ne2-ce1-he1: ?thet @exp @dif
*

!cg-nd1-hd1
set exp 126.88
quick @cg @nd1 @hd1
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*cg-nd1-hd1: ?thet @exp @dif
*

!ce1-nd1-hd1
set exp 126.22
quick @ce1 @nd1 @hd1
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ce1-nd1-hd1: ?thet @exp @dif
*

!cd2-cg-hg
set exp 132.60
quick @cd2 @cg @hg
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*cd2-cg-hg: ?thet @exp @dif
*

!nd1-cg-hg
set exp 121.92
quick @nd1 @cg @hg
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*nd1-cg-hg: ?thet @exp @dif
*

!ne2-cd2-hd2
set exp 121.36
quick @ne2 @cd2 @hd2
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*ne2-cd2-hd2: ?thet @exp @dif
*

!cg-cd2-hd2
set exp 127.95
quick @cg @cd2 @hd2
calc dif = ?thet - @exp
calc adif = @adif + abs ( @dif )
calc adif2 = @adif2 + (@dif * @dif)

calc count = @count + 1

write title unit 90
*cg-cd2-hd2: ?thet @exp @dif
*

!calculate statistics
! average difference
calc avediff = @adif / @count
calc avediff2 = @adif2 / @count
calc stdev = sqrt ( ( @avediff2 - ( @avediff * @avediff ) ) / (@count - 1) )

write title unit 90
* angle difference statistics
* average difference: @avediff
* standard deviation: @stdev
* ____________________________
*


stop