{+ file: bdomain_twin.inp +} {+ directory: xtal_twin +} {+ description: Grouped, unrestrained B-factor refinement of large domains for data with hemihedral twinning +} {+ authors: Axel T. Brunger, and Paul D. Adams +} {+ copyright: Yale University +} {+ reference: T.O. Yeates, Detecting and Overcoming Crystal Twinning, Methods in Enzymology 276, 344-358 (1997) +} {+ 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 -} {========================== IMPORTANT NOTE ===========================} {* This input file requires that the test set has been generated with the CNS task file make_cv_twin.inp. This ensures that twin related reflections are part of the same set. *} {- begin block parameter definition -} define( {======================= molecular structure =========================} {* molecular topology file *} {===>} structure_infile="porin.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="porin.pdb"; {====================== crystallographic data ========================} {* space group *} {* use International Table conventions with subscripts substituted by parenthesis *} {===>} sg="R3"; {* unit cell parameters in Angstroms and degrees *} {+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +} {===>} a=104.400; {===>} b=104.400; {===>} c=124.250; {===>} alpha=90; {===>} beta=90; {===>} gamma=120; {* 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="porin.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.25; {* 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=""; {============ initial 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; {======================= twinning parameters =========================} {* twinning operation *} {===>} twin_oper="h,-h-k,-l"; {* twinning fraction *} {===>} twin_frac=0.304; {========================== atom selection ===========================} {* select atoms to be included in refinement *} {===>} atom_select=(known and not hydrogen); {* atom selections which will be refined as groups *} {* any unselected atoms will remain unchanged note: the selections must be non-overlapping *} {* select atoms in group 1, use the selection all if everything is to be a group *} {===>} bd_group_1=(all); {* select atoms in group 2, use the selection none if not required *} {===>} bd_group_2=(none); {* select atoms in group 3, use the selection none if not required *} {===>} bd_group_3=(none); {* select atoms in group 4, use the selection none if not required *} {===>} bd_group_4=(none); {* select atoms in group 5, use the selection none if not required *} {===>} bd_group_5=(none); {* select atoms in group 6, use the selection none if not required *} {===>} bd_group_6=(none); {* select atoms in group 7, use the selection none if not required *} {===>} bd_group_7=(none); {* select atoms in group 8, use the selection none if not required *} {===>} bd_group_8=(none); {* select atoms in group 9, use the selection none if not required *} {===>} bd_group_9=(none); {* select atoms in group 10, use the selection none if not required *} {===>} bd_group_10=(none); {================= B-factor minimization parameters ==================} {* reset all atomic B factors to this number if positive *} {===>} reset_b=-1; {* number of steps of group B-factor minimization *} {===>} bfactor_nstep=30; {* B-factor limits *} {+ table: rows=1 "B-factor" cols=2 "minimum" "maximum" +} {===>} bmin=1; {===>} bmax=200; {* refinement target *} {+ list: twin_lsq: least squares residual for hemihedral twinning +} {+ choice: "twin_lsq" +} {===>} reftarget="twin_lsq"; {* 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 coordinate file *} {===>} coordinate_outfile="bdomain_twin.pdb"; {* format output coordinates for use in o *} {* if false then the default CNS output coordinate format will be used *} {+ choice: true false +} {===>} pdb_o_format=true; {===========================================================================} { 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 set echo=off end if ( &twin_frac > 0.5 ) then display Error: twinning fraction must be less than or equal 0.5 abort end if set echo=on end {- 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 ( $obs_low <= amplitude(remap[&STRIP%twin_oper]($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 8/01/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 if ( &reset_b > 0 ) then do (b=&reset_b) (&atom_select) end if xray predict mode=reciprocal to=fcalc selection=(ref_active=1) atomselection=( &atom_select ) end end {- BEGIN MODIFICATION 8/01/06 -} @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; 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; 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; ) {- 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_twin (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; twin_oper=&STRIP%twin_oper; twin_frac=&twin_frac; 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_twin (fobs=&STRIP%obs_f; fcalc=fcalc; fpart=fbulk; twin_oper=&STRIP%twin_oper; twin_frac=&twin_frac; 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_bd_group_$group = true ) then show sum(1) ( &bd_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 xray optimize group bmin=&bmin bmax=&bmax evaluate ($counter=1) while ( $counter <= $ngroup ) loop sele if ( $size_$counter > 0 ) then b=(&atom_select and (&bd_group_$counter)) end if evaluate ($counter=$counter+1) end loop sele nstep=&bfactor_nstep drop=1.0 ? end end xray predict mode=reciprocal to=fcalc selection=(ref_active=1) atomselection=( &atom_select ) end @@CNS_XTALMODULE:calculate_r_twin (fobs=&STRIP%obs_f; fcalc=fcalc; fpart=fbulk; twin_oper=&STRIP%twin_oper; twin_frac=&twin_frac; sel=(ref_active=1); sel_test=(tst_active=1); print=true; output=OUTPUT; r=$full_r; test_r=$full_test_r;) end set display=&coordinate_outfile end display REMARK coordinates from twinned data domain B-refinement display REMARK twinning operator= &STRIP%twin_oper twinning fraction= &twin_frac display REMARK refinement resolution: &low_res - &high_res A if ( $total_test > 0 ) then display REMARK starting twinned r= $start_r[f6.4] twinned free_r= $start_test_r[f6.4] display REMARK final twinned r= $full_r[f6.4] twinned free_r= $full_test_r[f6.4] else display REMARK starting twinned r= $start_r[f6.4] display REMARK final twinned r= $full_r[f6.4] end if xray wa=? end evaluate ($wa_print=$result) display REMARK wa= $wa_print display REMARK target= &STRIP%reftarget steps= &bfactor_nstep display REMARK sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma if ( &BLANK%parameter_infile_1 = false ) then display REMARK parameter file 1 : &STRIP%parameter_infile_1 end if if ( &BLANK%parameter_infile_2 = false ) then display REMARK parameter file 2 : &STRIP%parameter_infile_2 end if if ( &BLANK%parameter_infile_3 = false ) then display REMARK parameter file 3 : &STRIP%parameter_infile_3 end if if ( &BLANK%parameter_infile_4 = false ) then display REMARK parameter file 4 : &STRIP%parameter_infile_4 end if if ( &BLANK%parameter_infile_5 = false ) then display REMARK parameter file 5 : &STRIP%parameter_infile_5 end if display REMARK molecular structure file: &STRIP%structure_infile display REMARK input coordinates: &STRIP%coordinate_infile if ( &BLANK%anom_library = false ) then display REMARK anomalous f' f'' library: &STRIP%anom_library end if if ( &BLANK%reflection_infile_1 = false ) then display REMARK reflection file= &STRIP%reflection_infile_1 end if if ( &BLANK%reflection_infile_2 = false ) then display REMARK reflection file= &STRIP%reflection_infile_2 end if if ( &BLANK%reflection_infile_3 = false ) then display REMARK reflection file= &STRIP%reflection_infile_3 end if if ( &BLANK%ncs_infile = false ) then display REMARK ncs= &STRIP%ncs_type ncs file= &STRIP%ncs_infile else display REMARK ncs= none end if if ( &bscale # "no" ) then if ( $low_b_flag = true ) then display REMARK warning: B-correction gave atomic B-values less than zero display REMARK they have been reset to zero end if end if ! ! Begin modification (8/01/06) if ( &bscale = "anisotropic" ) then display REMARK Anisotropic B-factor tensor Ucart of atomic model without isotropic component : display REMARK B11=$Baniso_11[f8.3] B22=$Baniso_22[f8.3] B33=$Baniso_33[f8.3] display REMARK B12=$Baniso_12[f8.3] B13=$Baniso_13[f8.3] B23=$Baniso_23[f8.3] display REMARK Isotropic component added to coordinate array B: $Biso_model[f8.3] elseif ( &bscale = "isotropic" ) then display REMARK B-factor applied to coordinate array B: $Biso_model[f8.3] else display REMARK initial B-factor correction: none end if if ( &bulk_sol = true ) then display REMARK bulk solvent: probe radius=$solrad_best, shrink value=$solrad_best display REMARK bulk solvent: density level= $sol_k_ref e/A^3, B-factor= $sol_b_ref A^2 else display REMARK bulk solvent: false end if {- END MODIFICATION -} if ( &obs_type = "intensity" ) then display REMARK reflections with Iobs/sigma_I < &sigma_cut rejected display REMARK reflections with Iobs > &obs_rms * rms(Iobs) rejected display REMARK reflections with Iobs[&STRIP%twin_oper] = 0 rejected else display REMARK reflections with |Fobs|/sigma_F < &sigma_cut rejected display REMARK reflections with |Fobs| > &obs_rms * rms(Fobs) rejected display REMARK reflections with |Fobs|[&STRIP%twin_oper] = 0 rejected end if xray anomalous=? end if ( $result = true ) then display REMARK anomalous diffraction data was input end if {- MODIFIED 8/01/06 -} display REMARK fft gridding factor = $fft_grid, B factor offset = $fft_b_add A^2, Elimit = $fft_elim {- END MODIFICATION -} display REMARK theoretical total number of refl. in resol. range: $total_theor[I6] ( 100.0 % ) display REMARK number of unobserved reflections (no entry or |F|=0): $unobserved[I6] ( $per_unobs[f5.1] % ) display REMARK number of reflections rejected: $rejected[I6] ( $per_reject[f5.1] % ) display REMARK total number of reflections used: $total_used[I6] ( $per_used[f5.1] % ) display REMARK number of reflections in working set: $total_work[I6] ( $per_work[f5.1] % ) display REMARK number of reflections in test set: $total_test[I6] ( $per_test[f5.1] % ) remark @CNS_XTALMODULE:write_pdb (pdb_o_format=&pdb_o_format; coordinate_outfile=&coordinate_outfile; sgparam=$sgparam;) stop