{+ file: fp_fdp_group.inp +} {+ directory: xtal_refine +} {+ description: Grouped, unrestrained refinement of f' and f'' for selected atoms +} {+ comment: This is intended to derive f' and f'' values for anomalous scatters which are part of an atomic model. The anomalous library produced can then be used in refinement of the model against anomalous data. For example: to determine the f' and f'' values for seleniums at the anomalous peak wavelength, followed by refinement of the model, including seleniums, against data at that peak +} {+ authors: Axel T. Brunger, and Paul D. Adams +} {+ copyright: Yale University +} {+ reference: A.T. Brunger, The Free R Value: a Novel Statistical Quantity for Assessing the Accuracy of Crystal Structures, Nature 355, 472-474 (1992) +} {- Guidelines for using this file: - all strings must be quoted by double-quotes - logical variables (true/false) are not quoted - do not remove any evaluate statements from the file -} {- begin block parameter definition -} define( {======================= molecular structure =========================} {* molecular topology file *} {===>} structure_infile="amy.mtf"; {* parameter files *} {===>} parameter_infile_1="CNS_TOPPAR:protein_rep.param"; {===>} parameter_infile_2=""; {===>} parameter_infile_3=""; {===>} parameter_infile_4=""; {===>} parameter_infile_5=""; {* coordinate file *} {===>} coordinate_infile="amy.pdb"; {====================== crystallographic data ========================} {* space group *} {* use International Table conventions with subscripts substituted by parenthesis *} {===>} sg="P2(1)2(1)2(1)"; {* unit cell parameters in Angstroms and degrees *} {+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +} {===>} a=61.76; {===>} b=40.73; {===>} c=26.74; {===>} alpha=90; {===>} beta=90; {===>} gamma=90; {* anomalous f' f'' library file *} {* If a file is not specified, no anomalous contribution will be included *} {+ choice: "CNS_XRAYLIB:anom_cu.lib" "CNS_XRAYLIB:anom_mo.lib" "" user_file +} {===>} anom_library=""; {* reflection files *} {* specify non-anomalous reflection files before anomalous reflection files. *} {* files must contain unique array names otherwise errors will occur *} {===>} reflection_infile_1="amy.cv"; {===>} reflection_infile_2=""; {===>} reflection_infile_3=""; {* reciprocal space array containing observed amplitudes: required *} {===>} obs_f="fobs"; {* reciprocal space array containing sigma values for amplitudes: required *} {===>} obs_sigf="sigma"; {* reciprocal space array containing test set for cross-validation: required *} {* cross-validation should always be used, with the possible exception of a final round of refinement including all data *} {* cross-validation is always required for the maximum likelihood targets *} {===>} test_set="test"; {* number for selection of test reflections: required for cross-validation *} {* ie. reflections with the test set array equal to this number will be used for cross-validation, all other reflections form the working set *} {===>} test_flag=1; {* reciprocal space array containing weighting scheme for observed amplitudes: optional *} {* only used for the "residual" and "vector" targets - this will default to a constant value of 1 if array is not present *} {===>} obs_w=""; {* reciprocal space array containing observed intensities: optional *} {* required for the "mli" target *} {===>} obs_i=""; {* reciprocal space array containing sigma values for intensities: optional *} {* required for the "mli" target *} {===>} obs_sigi=""; {* reciprocal space arrays with experimental phase probability distribution: optional *} {* Hendrickson-Lattman coefficients A,B,C,D *} {* required for the "mlhl" target *} {+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +} {===>} obs_pa=""; {===>} obs_pb=""; {===>} obs_pc=""; {===>} obs_pd=""; {* complex reciprocal space array containing experimental phases: optional *} {* required for the "mixed" and "vector" targets *} {===>} obs_phase=""; {* reciprocal space array containing experimental figures of merit: optional *} {* required for the "mixed" target *} {===>} obs_fom=""; {* resolution limits to be used in refinement *} {* the full resolution range of observed data should be used in refinement. A bulk solvent correction should be applied to allow the use of low resolution terms. If no bulk solvent correction is applied, data must be truncated at a lower resolution limit of between 8 and 6 Angstrom. *} {+ table: rows=1 "resolution" cols=2 "lowest" "highest" +} {===>} low_res=500.0; {===>} high_res=2.0; {* apply rejection criteria to amplitudes or intensities *} {+ choice: "amplitude" "intensity" +} {===>} obs_type="amplitude"; {* Observed data cutoff criteria: applied to amplitudes or intensities *} {* reflections with magnitude(Obs)/sigma < cutoff are rejected. *} {===>} sigma_cut=0.0; {* rms outlier cutoff: applied to amplitudes or intensities *} {* reflections with magnitude(Obs) > cutoff*rms(Obs) will be rejected *} {===>} obs_rms=10000; {=================== non-crystallographic symmetry ===================} {* NCS-restraints/constraints file *} {* see auxiliary/ncs.def *} {===>} ncs_infile=""; {============ overall B-factor and bulk solvent corrections ==========} {* overall B-factor correction *} {+ choice: "no" "isotropic" "anisotropic" +} {===>} bscale="anisotropic"; {* bulk solvent correction *} {* a mask is required around the molecule(s). The region outside this mask is the solvent region *} {+ choice: true false +} {===>} bulk_sol=true; {* bulk solvent mask file *} {* mask will be read from O type mask file if a name is given otherwise calculated from coordinates of selected atoms *} {===>} bulk_mask_infile=""; {* automatic bulk solvent parameter search *} {+ choice: true false +} {===>} sol_auto=true; {* optional file with a listing of the results of the automatic bulk solvent grid search *} {===>} sol_output=""; {* fixed solvent mask parameters if the automatic option is not used *} {+ table: rows=1 "bulk solvent" cols=2 "probe radius (A)" "shrink radius (A)" +} {===>} sol_rad=1.0; {===>} sol_shrink=1.0; {* fixed solvent parameters if the automatic option is not used *} {+ table: rows=1 "bulk solvent" cols=2 "e-density level (e/A^3)" "B-factor (A^2)" +} {===>} sol_k=-1; {===>} sol_b=-1; {========================== atom selection ===========================} {* select atoms to be included in refinement *} {===>} atom_select=(known and not hydrogen); {* atom selections which will have f'/f'' refined as groups *} {* any unselected atoms will remain unchanged note: the selections must be non-overlapping *} {* selection should normally be on the basis of atom type: (chemical Se) - all seleniums (chemical Zn*) - all zincs *} {* select atoms in group 1, use the selection none if not required *} {===>} anom_group_1=(chemical S); {* select atoms in group 2, use the selection none if not required *} {===>} anom_group_2=(none); {* select atoms in group 3, use the selection none if not required *} {===>} anom_group_3=(none); {* select atoms in group 4, use the selection none if not required *} {===>} anom_group_4=(none); {=================== f'/f'' minimization parameters ==================} {* number of steps of group f'/f''minimization *} {===>} anom_nstep=30; {* minimum and maximum allowed f' *} {===>} fp_min=-100; {===>} fp_max=100; {* minimum and maximum allowed f'' *} {===>} fdp_min=-100; {===>} fdp_max=100; {* refinement target *} {+ list: mlf: maximum likelihood target using amplitudes mli: maximum likelihood target using intensities mlhl: maximum likelihood target using amplitudes and phase probability distribution residual: standard crystallographic residual vector: vector residual mixed: (1-fom)*residual + fom*vector e2e2: correlation coefficient using normalized E^2 e1e1: correlation coefficient using normalized E f2f2: correlation coefficient using F^2 f1f1: correlation coefficient using F +} {+ choice: "mlf" "mli" "mlhl" "residual" "vector" "mixed" "e2e2" "e1e1" "f2f2" "f1f1" +} {===>} reftarget="mlf"; {* number of bins for refinement target *} {* this will be determined automatically if a negative value is given otherwise the specified number of bins will be used *} {===>} target_bins=-1; {* memory allocation for FFT calculation *} {* this will be determined automatically if a negative value is given otherwise the specified number of words will be allocated *} {===>} fft_memory=-1; {=========================== output files ============================} {* output anomalous scattering factor library file *} {===>} anomlib_outfile="fp_fdp_group.lib"; {===========================================================================} { things below this line do not normally need to be changed } {===========================================================================} ) {- end block parameter definition -} checkversion 1.2 evaluate ($log_level=quiet) structure @&structure_infile end coordinates @&coordinate_infile parameter if ( &BLANK%parameter_infile_1 = false ) then @@¶meter_infile_1 end if if ( &BLANK%parameter_infile_2 = false ) then @@¶meter_infile_2 end if if ( &BLANK%parameter_infile_3 = false ) then @@¶meter_infile_3 end if if ( &BLANK%parameter_infile_4 = false ) then @@¶meter_infile_4 end if if ( &BLANK%parameter_infile_5 = false ) then @@¶meter_infile_5 end if end xray @CNS_XTALLIB:spacegroup.lib (sg=&sg; sgparam=$sgparam;) a=&a b=&b c=&c alpha=&alpha beta=&beta gamma=&gamma @CNS_XRAYLIB:scatter.lib if ( &BLANK%reflection_infile_1 = false ) then reflection @@&reflection_infile_1 end end if if ( &BLANK%reflection_infile_2 = false ) then reflection @@&reflection_infile_2 end end if if ( &BLANK%reflection_infile_3 = false ) then reflection @@&reflection_infile_3 end end if end if ( &BLANK%anom_library = false ) then @@&anom_library else set echo=off end xray anomalous=? end if ( $result = true ) then display Warning: no anomalous library has been specified display no anomalous contribution will used in refinement end if set echo=on end end if {- copy define parameters of optional arrays into symbols so we can redefine them -} evaluate ($obs_i=&obs_i) evaluate ($obs_sigi=&obs_sigi) evaluate ($obs_w=&obs_w) xray @@CNS_XTALMODULE:checkrefinput ( reftarget=&reftarget; obs_f=&obs_f; obs_sigf=&obs_sigf; test_set=&test_set; obs_pa=&obs_pa; obs_pb=&obs_pb; obs_pc=&obs_pc; obs_pd=&obs_pd; obs_phase=&obs_phase; obs_fom=&obs_fom; obs_w=$obs_w; obs_i=$obs_i; obs_sigi=$obs_sigi; ) query name=fcalc domain=reciprocal end if ( $object_exist = false ) then declare name=fcalc domain=reciprocal type=complex end end if declare name=fbulk domain=reciprocal type=complex end do (fbulk=0) ( all ) binresolution &low_res &high_res mapresolution &high_res if ( &obs_type = "intensity" ) then if ( &BLANK%obs_i = true ) then display Error: observed intensity array is undefined display aborting script abort end if evaluate ($reject_obs=&obs_i) evaluate ($reject_sig=&obs_sigi) show min (amplitude(&STRIP%obs_i)) (all) evaluate ($obs_lower_limit=$result-0.1) else evaluate ($reject_obs=&obs_f) evaluate ($reject_sig=&obs_sigf) evaluate ($obs_lower_limit=0) end if declare name=ref_active domain=reciprocal type=integer end declare name=tst_active domain=reciprocal type=integer end do (ref_active=0) ( all ) do (ref_active=1) ( ( amplitude($STRIP%reject_obs) > $obs_lower_limit ) and ( &low_res >= d >= &high_res ) ) statistics overall completeness selection=( ref_active=1 ) end evaluate ($total_compl=$expression1) show sum(1) ( ref_active=1 ) evaluate ($total_read=$select) evaluate ($total_theor=int(1./$total_compl * $total_read)) show rms (amplitude($STRIP%reject_obs)) ( ref_active=1 ) evaluate ($obs_high=$result*&obs_rms) show min (amplitude($STRIP%reject_obs)) ( ref_active=1 ) evaluate ($obs_low=$result) do (ref_active=0) ( all ) do (ref_active=1) ( ( amplitude($STRIP%reject_obs) >= &sigma_cut*$STRIP%reject_sig ) and ( $STRIP%reject_sig # 0 ) and ( $obs_low <= amplitude($STRIP%reject_obs) <= $obs_high ) and ( &low_res >= d >= &high_res ) ) do (tst_active=0) (all) if ( &BLANK%test_set = false ) then do (tst_active=1) (ref_active=1 and &STRIP%test_set=&test_flag) end if show sum(1) ( ref_active=1 and tst_active=0 ) evaluate ($total_work=$select) show sum(1) ( ref_active=1 and tst_active=1 ) evaluate ($total_test=$select) evaluate ($total_used=$total_work+$total_test) evaluate ($unobserved=$total_theor-$total_read) evaluate ($rejected=$total_read-$total_used) evaluate ($per_unobs=100*($unobserved/$total_theor)) evaluate ($per_reject=100*($rejected/$total_theor)) evaluate ($per_used=100*($total_used/$total_theor)) evaluate ($per_work=100*($total_work/$total_theor)) evaluate ($per_test=100*($total_test/$total_theor)) associate fcalc ( &atom_select ) tselection=( ref_active=1 ) cvselection=( tst_active=1 ) method=FFT {- MODIFIED 2/15/06 -} end show min ( b ) ( &atom_select ) evaluate ($b_min=$result) @@CNS_XTALMODULE:fft_parameter_check ( d_min=&high_res; b_min=$b_min; grid=auto; fft_memory=&fft_memory; fft_grid=$fft_grid; fft_b_add=$fft_b_add; fft_elim=$fft_elim; ) xray {- END MODIFICATION -} tolerance=0.0 lookup=false end igroup interaction (&atom_select) (&atom_select) end flags exclude * include xref end if ( &BLANK%ncs_infile = false ) then inline @&ncs_infile end if xray predict mode=reciprocal to=fcalc selection=(ref_active=1) atomselection=( &atom_select ) end end {- BEGIN MODIFICATION -} @CNS_XTALMODULE:scale_and_solvent_grid_search ( bscale=&bscale; sel=( ref_active=1 ); sel_test=( tst_active=1 ); atom_select=( &atom_select ); bulk_sol=&bulk_sol; bulk_mask=&bulk_mask_infile; bulk_atoms=( &atom_select ); sol_auto=&sol_auto; sol_k=&sol_k; sol_b=&sol_b; sol_rad=&sol_rad; sol_shrink=&sol_shrink; fcalc=fcalc; obs_f=&STRIP%obs_f; obs_sigf=&STRIP%obs_sigf; obs_i=$STRIP%obs_i; obs_sigi=$STRIP%obs_sigi; fpart=fbulk; ! ! Begin modification (6/28/06) Baniso_11=$Baniso_11; Baniso_22=$Baniso_22; Baniso_33=$Baniso_33; Baniso_12=$Baniso_12; Baniso_13=$Baniso_13; Baniso_23=$Baniso_23; Biso=$Biso_model; ! End modification ! sol_k_best=$sol_k_ref; sol_b_best=$sol_b_ref; solrad_best=$solrad_best; shrink_best=$shrink_best; b=b; low_b_flag=$low_b_flag; sol_output=&sol_output; ) xray @@CNS_XTALMODULE:calculate_r ( fobs=&STRIP%obs_f; fcalc=fcalc; fpart=fbulk; sel=( ref_active=1 ); sel_test=( tst_active=1 ); print=true; output=OUTPUT; r=$start_r; test_r=$start_test_r;) end {- check the gridding again since the minimum B-factor may have changed -} show min ( b ) ( &atom_select ) evaluate ($b_min=$result) @@CNS_XTALMODULE:fft_parameter_check ( d_min=&high_res; b_min=$b_min; grid=auto; fft_memory=&fft_memory; fft_grid=$fft_grid; fft_b_add=$fft_b_add; fft_elim=$fft_elim; ) {- END MODIFICATION -} xray @@CNS_XTALMODULE:refinementtarget (target=&reftarget; sig_sigacv=0.07; mbins=&target_bins; fobs=&STRIP%obs_f; sigma=&STRIP%obs_sigf; weight=$STRIP%obs_w; iobs=$STRIP%obs_i; sigi=$STRIP%obs_sigi; test=tst_active; fcalc=fcalc; fpart=fbulk; pa=&STRIP%obs_pa; pb=&STRIP%obs_pb; pc=&STRIP%obs_pc; pd=&STRIP%obs_pd; phase=&STRIP%obs_phase; fom=&STRIP%obs_fom; sel=(ref_active=1); sel_test=(tst_active=1); statistics=true;) end @@CNS_XTALMODULE:one_term_wa (wa=$wa;) xray @@CNS_XTALMODULE:calculate_r (fobs=&STRIP%obs_f; fcalc=fcalc; fpart=fbulk; sel=(ref_active=1); sel_test=(tst_active=1); print=true; output=OUTPUT; r=$start_r; test_r=$start_test_r;) end evaluate ($ngroup=1) evaluate ($group=1) evaluate ($done=false) while ( $done = false ) loop group if ( &exist_anom_group_$group = true ) then show sum(1) ( &anom_group_$group ) evaluate ($size_$group=$result) evaluate ($ngroup=$ngroup+1) else evaluate ($done=true) evaluate ($ngroup=$ngroup-1) end if evaluate ($group=$group+1) end loop group evaluate ($counter=1) while ( $counter <= $ngroup ) loop sele if ( $size_$counter > 0 ) then do (scatter_fp=-5) (&anom_group_$counter) do (scatter_fdp=5) (&anom_group_$counter) end if evaluate ($counter=$counter+1) end loop sele xray optimize group fpmin=&fp_min fpmax=&fp_max fdpmin=&fdp_min fdpmax=&fdp_max evaluate ($counter=1) while ( $counter <= $ngroup ) loop sele if ( $size_$counter > 0 ) then fp=(&atom_select and (&anom_group_$counter)) fdp=(&atom_select and (&anom_group_$counter)) end if evaluate ($counter=$counter+1) end loop sele nstep=&anom_nstep drop=1.0 ? end end xray predict mode=reciprocal to=fcalc selection=(ref_active=1) atomselection=( &atom_select ) end @@CNS_XTALMODULE:calculate_r (fobs=&STRIP%obs_f; fcalc=fcalc; fpart=fbulk; sel=(ref_active=1); sel_test=(tst_active=1); print=true; output=OUTPUT; r=$full_r; test_r=$full_test_r;) end set display=&anomlib_outfile end display !=================================================================== display ! file: &anomlib_outfile automatically generated by fp_fdp_group.inp if ( $total_test > 0 ) then display ! starting r= $start_r[f6.4] free_r= $start_test_r[f6.4] display ! final r= $full_r[f6.4] free_r= $full_test_r[f6.4] else display ! starting r= $start_r[f6.4] display ! final r= $full_r[f6.4] end if display !=================================================================== display ! f'/f'' refined for the following groups of atoms: evaluate ($counter=1) while ( $counter <= $ngroup ) loop sele if ( $size_$counter > 0 ) then display ! --> &anom_group_$counter end if evaluate ($counter=$counter+1) end loop sele display evaluate ($counter=1) while ( $counter <= $ngroup ) loop sele if ( $size_$counter > 0 ) then show ave(scatter_fp) (&atom_select and &anom_group_$counter) display do (scatter_fp=$result) &anom_group_$counter show ave(scatter_fdp) (&atom_select and &anom_group_$counter) display do (scatter_fdp=$result) &anom_group_$counter display end if evaluate ($counter=$counter+1) end loop sele set display=OUTPUT end stop