! refine.inp -- Clean up the output of dgsa.inp using simulated annealing,
!     including Van der Waals and electrostatic forces
! Dave Schweisguth <dcs@proton.chem.yale.edu>, 21 Jul 1996
! Derived from nmr/refine_gentle.inp

evaluate ($t=300)		! Temperature for constant-temperature MD
evaluate ($steps=40000)    	! Number of steps of constant-temperature MD
evaluate ($timestep=0.0005)	! Time of each MD step

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
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

coordinates @dgsa.pdb		! END in PDB file

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

! Intial minimization.

minimize powell nstep=200 nprint=25 end

! Constant temperature SA

vector do (fbeta=100) (all)	! Coupling to heat bath.
dynamics verlet
    iasvel=maxwell
    firsttemperature=$t
    tbath=$t
    tcoup=true
    nstep=$steps
    timestep=$timestep
    nprint=250
    iprfrq=2500
    trajectory=refine.tra
    nsavc=200			! Save trajectory to file.
end
close refine.tra end

! Compute dynamics average over second half of trajectory.

evaluate ($steps2=$steps/2)
dynamics analysis average
    input=refine.tra
    begin=$steps2
    skip=200
    stop=$steps
end
close refine.tra end		! disposition=delete doesn't work

! Final minimization.

minimize powell nstep=200 nprint=25 end

! Analyze and write out the final 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

write coordinates output=refine.pdb end
    
stop