! average.inp -- Computes the average structure, RMSDs from the average
! for each original structure and the rms RMSD.
! Dave Schweisguth <dcs@proton.chem.yale.edu>, 1 Sep 1996
! Derived from nmr/average.inp

set seed=@xplor.seed end	! Use 'xplor -s'

set echo=off message=off end	! Normal use
!set echo=on message=all end	! Auxiliary file debugging

structure @rsf.inp end		! Read the structure file.
parameter @XPLORLIB:toppar/dna-rna-allatom.par end

noe nres=1000 @noes.stem.inp @noes.loop.inp @base-pairs.inp end
restraints dihedral nassign=300 @dihedrals.inp end
evaluate($planarweight=50)	! AAFS' choice
restraints planar @planars.inp end

set echo=on message=all end

noe				! Parameters for NOE effective energy term.
    ceiling=1000                      
    averaging  * cent
    potential  * square
    sqconstant * 1
    sqexponent * 2
    scale      * 100		! Constant NOE scale throughout the protocol.
end

parameter
    nbonds
	ctonnb=9.5
	ctofnb=10.5
	cutnb=11.5
	tolerance=0.5
	rdie
	switch
	vswitch
    end
end

restraints dihedral scale=800 end

flags exclude * include bond angl dihe impr vdw elec noe cdih plan end

! Compute the mean structure.

evaluate ($nfile=9)
vector do (store1=0) (all)
vector do (store2=0) (all)
vector do (store3=0) (all)
vector do (store4=0) (all)

evaluate ($count=0)
while ($count < $nfile) loop main
    evaluate ($count=$count+1)
    evaluate ($filename="average."+encode($count)+".pdb")
    coordinates initialize end
    coordinates @@$filename
    if ($count=1) then
	coordinates copy end	! Store first structure in comparison set.
    end if
    coordinates fit end
    vector do (store1=store1+x) (all)
    vector do (store2=store2+y) (all)
    vector do (store3=store3+z) (all)
    vector do (store4=store4+x*x+y*y+z*z) (all)
end loop main

vector do (x=store1 / $nfile) (all)
vector do (y=store2 / $nfile) (all)
vector do (z=store3 / $nfile) (all)
vector do (b=sqrt(max(0,store4/$nfile-(x**2+y**2+z**2)))) (all)

! Minimize the averaged structure

minimize powell nstep=200 nprint=25 end

! Analyze and write out the averaged structure

print threshold=0.05 bonds          
evaluate ($rms_bond=$result)
evaluate ($v_bond=$violations)
print threshold=5 angles
evaluate ($rms_angl=$result)
evaluate ($v_angl=$violations)
print threshold=15 dihedrals
evaluate ($rms_dihe=$result)
evaluate ($v_dihe=$violations)
print threshold=5 impropers
evaluate ($rms_impr=$result)
evaluate ($v_impr=$violations)
print threshold=0.2 noe 
evaluate ($rms_noe=$result)
evaluate ($v_noe=$violations)
print threshold=5 cdih
evaluate ($rms_cdih=$result)
evaluate ($v_cdih=$violations)

set echo=off message=off end

display Energy:      bond $bond, angle $angl, dihedral $dihe,
display              improper $impr, NOE $noe, c-dihedral $cdih,
display              planar $plan, VdW $vdw, elec $elec, total $ener
display RMSD:        bond $rms_bond, angle $rms_angl, dihedral $rms_dihe,
display              improper $rms_impr, NOE $rms_noe, c-dihedral $rms_cdih
display Violations:  bond $v_bond, angle $v_angl, dihedral $v_dihe,
display              improper $v_impr, NOE $v_noe, c-dihedral $v_cdih

set echo=on message=all end

! b array (last column) is the rms difference from the mean

write coordinates output=average.pdb end

! Compute an overall rms difference.

coordinates copy end
evaluate ($rmsd_rmsd=0)
evaluate ($rmsd_rmsd_stem=0)
evaluate ($rmsd_rmsd_loop=0)

evaluate ($count=0)
while ($count < $nfile) loop main
    evaluate ($count=$count+1)
    evaluate ($filename="average."+encode($count)+".pdb")
    coordinates initialize end
    coordinates @@$filename
    coordinates fit end
    coordinates rms end
    evaluate ($rmsd_rmsd=$rmsd_rmsd + $result**2)
    coordinates fit selection=(residue 1:7 or residue 15:21) end
    coordinates rms selection=(residue 1:7 or residue 15:21) end
    evaluate ($rmsd_rmsd_stem=$rmsd_rmsd_stem + $result**2)
    coordinates fit selection=(residue 8:14) end
    coordinates rms selection=(residue 8:14) end
    evaluate ($rmsd_rmsd_loop=$rmsd_rmsd_loop + $result**2)
end loop main

! RMS RMSD to the average structure

evaluate ($rmsd_rmsd=sqrt($rmsd_rmsd / $nfile))
evaluate ($rmsd_rmsd_stem=sqrt($rmsd_rmsd_stem / $nfile))
evaluate ($rmsd_rmsd_loop=sqrt($rmsd_rmsd_loop / $nfile))

stop