* general drug molecule, minimize and write
* input for Gaussian QM package
*

!set identifier
set id1 drug_1
set toppar $CHARMMDATA

set top @toppar/top_all27_prot_lipid.rtf
set par @toppar/par_all27_prot_lipid.prm
set toppar top_mmtsb_example.str

! tolgrd setting: RMS force tolerance to exit minimizer
set a 0.0001

! read the topology
open unit 9 read form name @top
read rtf card unit 9
close unit 9

! read the parameters
open unit 9 read form name @par
read para card unit 9
close unit 9

!stream toppar file with additional rtf/param information
stream @toppar

read sequence cards
* cmpd model 2
*
1
mod2

generate a first none last none setup warn

!generate preliminary IC table
!ic generate
!stop

ic param
ic seed 1 cg 1 cd1 1 ne1
ic build
ic print

coor print

!write initial coordinates
open unit 26 write form name @id1_gen.pdb
write coor pdb unit 26

!minimization with infinite nonbond cutoff
update cutnb 99.0 ctofnb 98.0 ctonnb 97.0 switch vswitch
mini abnr nstep 200 nprint 10
mini nraph nstep 50

!write minmized coordinates
open unit 26 write form name @id1_min.pdb
write coor pdb unit 26

ic fill
ic print

!Create input file for Gaussian QM program. Includes both
!energy minimization and frequency calculation.  Note the
!total charge is assummed to be zero
!
!mp2 calculations on the full drug molecule will be difficult,
!therefore HF may be more appropriate
!

open unit 30 write form name gauss/@id1_opt_freq_mp2.com
prnlev 3
mmqm unit 30
%chk=@id1_opt_freq_mp2.chk
%nproc=2
#MP2/6-31G* Opt pop=(mk,dipole) density=current

@id1 optimized at MP2/6-31G*

0 1
END
GAUSSIAN_BASIS


--link1--
%nprocs=2
%chk=@id1_opt_freq_mp2.chk
#MP2/6-31G* geom=allcheck freq IOP(7/33=1)

@id1 frequency

0 1


END !end write gaussian input file

stop