!
! 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_KSOL = 0.65)
!
! We just set the value of KSOL to be used in this run.
! You could write another script to loop over different values of
! $KSOL, or just edit this one and run it again with a new $TRIAL_KSOL.
!
! 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 ($KSOL = $TRIAL_KSOL)

	eval ($BSOL = 20.)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_B20 = $RESULT)
	end
	vector do (B = BCOMP) (all)

	eval ($BSOL = 40.)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_B40 = $RESULT)
	end
	vector do (B = BCOMP) (all)

	eval ($BSOL = 60.)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_B60 = $RESULT)
	end
	vector do (B = BCOMP) (all)

	eval ($BSOL = 80.)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_B80 = $RESULT)
	end
	vector do (B = BCOMP) (all)

	eval ($BSOL = 100.)
	@bulksol.inp
	xrefine
	    print R-factor
	    eval ($R_B100 = $RESULT)
	end
	vector do (B = BCOMP) (all)

set display=bsolref.dat end
display # Search for optimal BSOL with KSOL = $TRIAL_KSOL
display # BSOL   R-value     (orig R = $R_orig)
display 20. $R_B20
display 40. $R_B40
display 60. $R_B60
display 80. $R_B80
display 100. $R_B100

stop