remarks file geomanal/ramachandran.inp -- Make a Ramachandran (phi-psi) plot
remarks      (this file is specific for proteins)

                                                     {*Read structure file.*}
{===>} structure @../generate/generate.psf end 
 
                                                        {*Read coordinates.*}
{===>} coordinates @../xtalrefine/slowcool.pdb    

parameter 
   {===>} @TOPPAR:parhcsdx.pro          {*Read empirical energy parameters.*}
end

{===>} vector ident (store9) ( name ca )     {*This selects all residues to*}
                                             {*be analyzed.   Change it    *}
                                             {*to analyze portions, e.g.,  *}
                                             {*( name ca and segid "A" ).  *}
                                           {*"tag" assigns one unique atom *}
                                           {*per residue (see <selection>).*}
                 
{===>} set display=ramachandran.list end    {*Write the phi psi info to the*}
                                            {*specified file.              *} 

set echo=off end                          {*Turn off echo to reduce output.*}
set message=off end                       {*Turn off warning messages.     *}


evaluate ($pre_atom_id=0)                {*Initialize atom ids for residue *}
evaluate ($pre_pre_atom_id=0)            {*n-1 and n-2.                    *}

for $atom_id in id ( name ca ) loop main          {*Loop over all ca atoms.*}
    if ( $pre_pre_atom_id > 0 ) then


        pick  dihedral                                          {*Pick phi.*}
            ( byresidue ( id $pre_pre_atom_id ) and name C  )
            ( byresidue ( id $pre_atom_id ) and name N  )
            ( byresidue ( id $pre_atom_id ) and name CA  )
            ( byresidue ( id $pre_atom_id ) and name C  )
             geometry
        evaluate ($phi=$result)

        pick  dihedral                                          {*Pick psi.*}
            ( byresidue ( id $pre_atom_id ) and name N  )
            ( byresidue ( id $pre_atom_id ) and name CA  )
            ( byresidue ( id $pre_atom_id ) and name C  )
            ( byresidue ( id $atom_id ) and name N  )
             geometry
        evaluate ($psi=$result)

        vector show element ( resname ) ( byresidue ( id $pre_atom_id ) 
                                        and name ca ) 
        evaluate ( $resname=$result ) 

                                           {*Depending on the residue name,*}
                                           {* write a "1" or a "0" .       *}
        if ($resname = "GLY" ) then
           display    0  $phi $psi 
        else
           display    1  $phi $psi
        end if

    end if
    evaluate ($pre_pre_atom_id=$pre_atom_id)           {*Reassign atom ids.*}
    evaluate ($pre_atom_id=$atom_id) 
end loop main


stop