* This file is interface.inp.
* Usage, charmm <interface.inp >interface.out.
* Input files, top_all27_prot_na.rtf, par_all27_prot_na.prm, subunit1.pdb, subunit2.pdb from
*  fixpdb.awk.
* Reads proteins, adds H atoms, adds missing side chain atoms, determines atoms of subunit 1
*  that line within 4 Angstroms of any atom of subunit 2
* 

! Open and read topology and parameter files
open read card unit 20 name "top_all27_prot_na.rtf" 
read rtf card unit 20
close unit 20

open read card unit 20 name "par_all27_prot_na.prm" 
read parameter card unit 20
close unit 20

! Read in the subunits
open read unit 21 card name "subunit1.pdb" 
read sequence pdb unit 21

generate pro1 setup 
rewind unit 21

open read unit 22 card name "subunit2.pdb"
read sequence pdb unit 22

generate pro2 setup
rewind unit 22

read coordinate pdb offset -6 unit 21
close unit 21

read coordinate pdb append offset -6 unit 22
close unit 22

! Build missing atoms
ic fill preserve
ic parameter
ic build
hbuild

open write unit 30 card name prot.pdb
write coordinates pdb select segid pro2 .and. (segid pro1 .around. 4) show end unit 30
* Coordinates of all atoms in protein, pdb format
*

stop 

! Write outputs
open write unit 30 card name prot.pdb
write coordinates pdb select all end unit 30
* Coordinates of all atoms in protein, pdb format
*

open write unit 31 card name prot.crd
write coordinates card unit 31
* Coordinates of all atoms in protein, crd format 
*

open write unit 32 card name prot.psf
write psf unit 32 card
* psf of the protein
*

coor stat select all end
system -
"echo Xmin=?XMIN Xmax=?XMAX $'\n'YMin=?YMIN Ymax=?YMAX $'\n'Zmin=?ZMIN Zmax=?ZMAX  > /dev/tty" 



stop