*Generate psf and point mutation for a trajectory with p coordinates and i residues.
*Assumes a generated trajectory present.
*Coordinates must be temporarily stored in charmm-format before removing
*and rebuilding missing atoms to an optional residue (i.e. valine to alanine or glycine)
*This example is based on simulations done on the glucocorticoid receptor dna-binding domain
*(Int. J. Quant. Chem., 2001, 83, 230 -, PROTEINS, 2002, 49, 24-)
*Johan Bredenberg - 2003
*
nobo

read rtf card name $CHM_HOME/toppar/top_all22_prot.inp
read para card name $CHM_HOME/toppar/par_all22_prot.inp
!first declare wt psf

read sequ card
*title:amino_acids
*
71
cyc leu val cyc ser asp glu ala ser gly -
cys hsp tyr gly val leu thr cyc gly ser -
cyc lys val phe phe lys arg ala val glu -
gly gln hsp asn tyr leu cyc ala gly arg -
asn asp cyc ile ile asp lys ile arg arg -
lys asn cyc pro ala cyc arg tyr arg lys -
cys leu gln ala gly met asn leu glu ala -
arg
generate dbd1 setup
read sequ zn 2
gene ion2 setup

scalar charge stat sele type ca end ! use ca since they're always there.
!number of residues to be mutated, those already being present i.e. ala to ala
!will not have to be mutated.
!in this example, however, all residues are included.


set totres ?nsel
set k = 1

!read in the wild-type generated trajectory

open unit 34 read unform name tmp.trj
trajec query unit 34
set n ?nfile
trajec firstu 34 nuni 1 begin 1

set p = 1 !counter for extracting the coordinates
!might be a good idea to implement a code for mutating on-the-fly

label dumpcoord
trajec read
coor stat
coor orie norot
coor stat
write coor card name @p.crd !dump physical charmm coordinates
incr p by 1
if p le @n goto dumpcoord
dele atom sele all end !remove everything

label INIT
set p = 1 ! re-initiate coord frame reading
set i @k !residue number

set @i = ala

!Make a psf for the mutation of residue "i" to alanine

read sequ card
*AMINO-ACIDS
*
71
@1 -
@2 -
@3 -
@4 -
@5 -
@6 -
@7 -
@8 -
@9 -
@10 -
@11 -
@12 -
@13 -
@14 -
@15 -
@16 -
@17 -
@18 -
@19 -
@20 -
@21 -
@22 -
@23 -
@24 -
@25 -
@26 -
@27 -
@28 -
@29 -
@30 -
@31 -
@32 -
@33 -
@34 -
@35 -
@36 -
@37 -
@38 -
@39 -
@40 -
@41 -
@42 -
@43 -
@44 -
@45 -
@46 -
@47 -
@48 -
@49 -
@50 -
@51 -
@52 -
@53 -
@54 -
@55 -
@56 -
@57 -
@58 -
@59 -
@60 -
@61 -
@62 -
@63 -
@64 -
@65 -
@66 -
@67 -
@68 -
@69 -
@70 -
@71
gene mut setup
read sequ zn 2
gene iom setup

!store the psf of the mutation

write psf card name @k.psf

!write a trajectory for the mutation by rebuilding the wild-type trajectory

open unit 35 write unform name @k.trj ! index number for mutation
trajec iwrite 35 nfile @n begin 1

label readc
read coor card name @p.crd ! coordinate number p from extracted wild-type trajectory
ic para
ic build
hbuild sele ( hydrogen .and. ( segid mut .and. ( resid @i ) ) ) end ! add beta-hydrogens
coor stat
coor orie mass
coor stat
trajec write
incr p by 1
if p le @n goto readc

rewind unit 34
close unit 35

!mutate next residue in sequence - or whatever -
!remove the temporary mutated psf

dele atom sele all end

incr k by 1
if k le @totres goto INIT

stop

!This should yield the modified trajectories which now can be used in !subsequent calulations.