!
! CYCLE.INP
!
! Ethan A Merritt
! This is my file for standard crystallographic refinement in XPLOR.
! Before each refinement edit the set of lines below to control the
! refinement behaviour.
! Also change the REMARK statements later in the script to properly
! identify your protein.
!

eval ($REFH = 0)	! non-zero to refine hydrogens first with all else fixed
eval ($BULK = 1)	! non-zero to apply bulk solvent correction
eval ($BREF = 0.0)	! non-zero to refine B values after adding GAUSS($BREF)
eval ($BANI = 1)	! non-zero to include an overall anisotropic B in FCALC
eval ($SHAKE= 0)	! non-zero to randomly perturb coords at start
eval ($GEOM = 1)	! monitor ideality of geometry

eval ($NCYCLES = 3)	! number of cycles of refinement

! This run:
!

eval	($CHAINID_IS_SEGID = 1)

 structure	@ generate.psf end
 coordinates	@ current.protein

!
! Add in discrete solvent model
!
  		@ solvent_model.inp
 coordinates	@ current.solvent

!
! Delete any atoms not yet in model
!
    delete	selection=(hydrogen and segid SOLV) end
    delete	selection=(not known) end

!
! Standard setup for refinement
! Note that the xprepare.inp file must set resolution limits, data cutoffs, etc.
! If you are using a bulk solvent model it must also set $KSOL and $BSOL.
!
	@xprepare.inp

xrefin
    update-fcalc
    print R-factor
end

{* Truncate B's before refining *}
if ($BREF # 0) then
    vector do (B = 45. + 0.25*(B-45.)) (attr B > 45.)
    vector do (B =  5. + 0.50*(B-5.))  (not hydrogen and attr B < 5.)
end if

{* Stuck in a rut - shake it up a tiny bit *}
if ($SHAKE # 0) then
    vector do (X = X + GAUSS(0.01)) (all)
    vector do (Y = Y + GAUSS(0.01)) (all)
    vector do (Z = Z + GAUSS(0.01)) (all)
end if

{* Check initial stereochemistry (protein only) *}
if ($GEOM #0) then
    vector ident (store9) (segid SOLV or hydrogen or (name CE and resname MET))
    constraint interactions (resid 1:103 and not recall 9)
			   =(resid 1:103 and not recall 9) end
    print threshold=1.00 bonds
    print threshold=180.0 angles
    constraint interactions (all)=(all) end
end if

{*====== Preliminary geometric optimization to move   ======*}
{*====== hydrogens and newly built atoms from conflict =====*}
if ($REFH # 0) then
    vector ident (store4) (name JUNK)

    constraint fix = (not hydrogen and not recall 4) end
    flags excl XREF end
  
    minimize powell
      nstep = 40 		! 40 cycles
      drop  = 10  	! expected initial drop in energy (not critical)
      end
  
    constraint fix = (name NONE) end
    flags incl XREF end
end if

{*====== Give the B's a good shake to get them out of ======*}
{*====== a rut;  later on we'll shake more gently     ======*}
if ($BREF # 0) then
   	vector do (B = B + GAUSS(2. * $BREF)) (not hydrogen)
end if


{*====== Main Loop over several rounds of refinement  ======*}
{*====== each time do positional first, then B's      ======*}
	evaluate ($cycle = 0)
	while ($cycle < $ncycles) loop MAIN   evaluate ($cycle = $cycle + 1)

!
! Start with a clean slate in Fpartial, Fcalc
!
    xrefine
	do complex (FPART = 0.0)
	update
    end
    
!
! Here comes our solvent correction
!
if ($BULK # 0) then
    eval ($KBSOLNAME = "")		! no extra output from here
    @bulksol.inp			! This loads correction into FPART
end if

!
! And here is an overall anisotropic B correction
!
if ($BANI #0) then
    @anisob.inp
end if

minimize powell
    nstep = 40 		! 40 cycles
    drop  = 10  	! expected initial drop in energy (not critical)
    ?
    end

if ($BREF # 0) then
    ! Shake up the B's to get out of a rut *}
   	vector do (B = B + GAUSS($BREF)) (not hydrogen)
    xrefin
  	optimize bfactors              	! invoke B-factor optimizer
	     nstep = 30                	! for lots of cycles
	     drop = 10.0               	! expected initial drop in energy
	     rweight = -1 		! set auto weight for restraints
	     bond-restraint = 2.0       ! target sigma for 1-2 B-factor pairs
	     angle-restraint = 4.0    	! target sigma for 1-3 (angle) pairs
	end
    end
end if

!xrefine
!	update
!	print R-factor
!end

{* Check protein stereochemistry at end of cycle *}
if ($GEOM #0) then
    vector ident (store9) (segid SOLV or hydrogen or (name CE and resname MET))
    constraint interactions (not segid SOLV and not recall 9)
			   =(not segid SOLV and not recall 9) end
    print threshold=1.00 bonds
    eval ($RMSBONDS = $RESULT)
    print threshold=180.0 angles
    eval ($RMSANGLES = $RESULT)
    constraint interactions (all)=(all) end
end if

if ($BULK # 0) then
    REMARK My protein name and vital statistics
    REMARK after positional + B refinement cycle $CYCLE   $DATE
    REMARK	R = $R     Rfree = $TEST_R
    REMARK non-ideality:	bonds = $RMSBONDS	angles = $RMSANGLES
    REMARK Bulk solvent model:		KSOL = $KSOL, BSOL = $BSOL
else
    REMARK My protein name and vital statistics
    REMARK after positional + B refinement cycle $CYCLE   $DATE
    REMARK	R = $R     Rfree = $TEST_R
    REMARK non-ideality:	bonds = $RMSBONDS	angles = $RMSANGLES
    REMARK Bulk solvent model:		KSOL = $KSOL, BSOL = $BSOL
end if

write coordinates 
	output=cycle.protein
	select= (not SEGID SOLV)
	end

write coordinates 
	output=cycle.solvent
	select= (SEGID SOLV)
	end

{*====== End of main loop ======*}
	end loop MAIN

if ($NCYCLES = 0) then
    write coordinates 
	output=cycle.protein
	select= (not SEGID SOLV)
    end
end if

xrefine
    update-fcalc
    print R-factor
    end

 xrefin
   wa = 10000.                          ! not at all clear whether this matters
   optimize bfactors                    ! invoke B-factor optimizer
      nstep = 1                         ! just to check targets
      rweight = -1                      ! auto-weighting of restraints
      drop = 10.0                       ! expected initial drop in energy
      bond-restraint = 1.0              ! target sigma for 1-2 B-factor pairs
      angle-restraint = 1.0             ! target sigma for 1-3 (angle) pairs
      end
   end

 stop