!
! Optimize choice of two parameters $KSOL and $BSOL for use in 
! bulk solvent correction using Babinet principal. 
! The correction to be applied is:
!
!		Fc = (1 - Ksol exp(-Bsol*s**2))Fc
!
! We will loop over several possible values of $BSOL, and dump the resulting
! Luzzati plots to a series of output files defined by the name
! result taken as indicating optimal choice of $BSOL for this value
! of $KSOL.
!
eval ($TRIAL_BSOL = 50.)
!
! We just set the value of BSOL to be used in this run.
! You could write another script to loop over different values of
! $BSOL, or just edit this one and run it again with a new $TRIAL_BSOL.
!
! Ethan A Merritt - 1995

eval	($CHAINID_IS_SEGID = 1)
eval	($KBSOLNAME = "")

structure	@ generate.psf end
coordinates	@ current.protein

{* Add in discrete solvent model *}
 		@ solvent_model.inp
coordinates	@ current.solvent

delete		selection=(not known) end

		@xprepare.inp

!
! Each run on the bulk solvent correction will fiddle with B factors
! so store the original ones to reset between trials
!
	vector do (BCOMP = B) (all)

!
! Document starting conditional
!
	eval ($KSOL = 0.0)
	eval ($BSOL = 0.0)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_orig = $RESULT)
	end
	vector do (B = BCOMP) (all)

!
! Try out bulk solvent model with several possible values of Bsol
!
	eval ($BSOL = $TRIAL_BSOL)

	eval ($KSOL = 0.5)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_K5 = $RESULT)
	end
	vector do (B = BCOMP) (all)

	eval ($KSOL = 0.6)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_K6 = $RESULT)
	end
	vector do (B = BCOMP) (all)

	eval ($KSOL = 0.7)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_K7 = $RESULT)
	end
	vector do (B = BCOMP) (all)

	eval ($KSOL = 0.8)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_K8 = $RESULT)
	end
	vector do (B = BCOMP) (all)

	eval ($KSOL = 0.9)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_K9 = $RESULT)
	end
	vector do (B = BCOMP) (all)

set display=ksolref.dat end
display # Search for optimal KSOL with BSOL = $TRIAL_BSOL
display # KSOL   R-value     (orig R = $R_orig)
display 0.5 $R_K5
display 0.6 $R_K6
display 0.7 $R_K7
display 0.8 $R_K8
display 0.9 $R_K9

stop