{+ file: cis_peptide.inp +}
{+ directory: general +}
{+ description: Identify cis-peptides in a protein structure +}
{+ comment: 
   writes out appropriate parameter file to maintain 
   cis-peptide geometry during refinement +}
{+ authors: Axel T. Brunger, and Paul D. Adams +}
{+ copyright: Yale University +}

{- 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 =========================}

{* structure file *}
{===>} structure_infile="fab2hfl.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="fab2hfl.pdb";

{========================== atom selection ===========================}

{* select atoms to be included in calculation *}
{===>} atom_select=(known and not hydrogen);

{=========================== output files ============================}

{* output parameter file to be included in subsequent refinement *}
{===>} list_outfile="cis_peptide.param";

{===========================================================================}
{        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
      @@&parameter_infile_1
   end if
   if ( &BLANK%parameter_infile_2 = false ) then
      @@&parameter_infile_2
   end if
   if ( &BLANK%parameter_infile_3 = false ) then
      @@&parameter_infile_3
   end if
   if ( &BLANK%parameter_infile_4 = false ) then
      @@&parameter_infile_4
   end if
   if ( &BLANK%parameter_infile_5 = false ) then
      @@&parameter_infile_5
   end if
 end

 if ( $log_level = verbose ) then
   set message=normal echo=on end
 else
   set message=off echo=off end
 end if

 evaluate ($first=true)
 evaluate ($written=false)

 for $id in id ( &atom_select and name ca ) loop main

   show (segid) (id $id)
   evaluate ($segid=$result)
   show (resid) (id $id)
   evaluate ($resid=$result)
   show (resname) (id $id)
   evaluate ($resname=$result)

   identity (store1) ( &atom_select and ( name c and bondedto 
            ( name n and resid $resid and segid $segid ) ) )
   if ( $select = 1 ) then
     show element (store1) (attribute store1 > 0)
     evaluate ($id_prev=$result)
     show (segid) (id $id_prev)
     evaluate ($segid_prev=$result)
     show (resid) (id $id_prev)
     evaluate ($resid_prev=$result)
     show (resname) (id $id_prev)
     evaluate ($resname_prev=$result)
 
     pick dihedral
       (name ca and segid $segid_prev and resid $resid_prev)
       (name  c and segid $segid_prev and resid $resid_prev)
       (name  n and segid $segid and resid $resid)
       (name ca and segid $segid and resid $resid)
       geometry
       
     evaluate ($dihedral=mod($result+360,360))

     if ( $dihedral > 180 ) then
       evaluate ($dihedral=$dihedral-360)
     end if

     evaluate ($absdihedral=abs($dihedral))

     set display=&list_outfile end

     if ( $absdihedral < 15 ) then
       if ( $first = true ) then
         display {===========================================================}
         display ! atom based parameter modifications to maintain 
         display !   cis-peptide geometry during refinement
         display ! automatically generated by script cis_peptide.inp
         display ! input structure file: &STRIP%structure_infile
         display ! input coordinate file: &STRIP%coordinate_infile
         display ! CNS version: $cns_version  date: $date  user: $name
         display {===========================================================}
         display 
         display set echo=off message=off end
         display 
         evaluate ($first=false)
         evaluate ($written=true)
       end if
       
       {- add double quotes to segid in order to allow blank segids -}
       evaluate ($segid="""+$segid+""") 
       evaluate ($segid_prev="""+$segid_prev+""")
       
       display {===========================================================}
       display !>>>> cis-peptide: segid=$segid resid=$resid resname=$resname  
       display !>>>>              current dihedral value= $dihedral[f8.3]
       display 
       display !to maintain this geometry use the 
       display !  following parameter statement during refinement:
       display
       display dihedral ( name ca and segid $segid_prev and resid $resid_prev )
       display          ( name  c and segid $segid_prev and resid $resid_prev )
       display          ( name  n and segid $segid and resid $resid )
       display          ( name ca and segid $segid and resid $resid ) 
       display            1250.0 1 180.0
       display
       display {===========================================================}
     elseif ( $absdihedral < 165 ) then
       display {===========================================================}
       display !>>>> distorted peptide: segid= $segid resid= $resid resname=$resname
       display !>>>>                    current dihedral value= $dihedral[f8.3]
       display 
       display !this may require correction before starting refinement
       display 
       display {===========================================================}
     end if

     set display=OUTPUT end
   else
     display couldn't find anything bonded to atom N
   end if

 end loop main

 if ( $written = true ) then
   set display=&list_outfile end
   display 
   display set echo=on message=on end
   display
   set display=OUTPUT end
 end if

 stop