remarks file geomanal/geomplot.inp -- Give conformational energy per residue
remarks Provide a printer output file as well as a list file 
remarks for plotting

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

  {===>} vector ident (store9) ( tag )    {*This selects all residues to   *}
                                          {*be analyzed.   Change it       *}
                                          {*to analyze portions,           *}
                                          {*e.g.,  ( tag and segid "A" ).  *}
                                          {*"tag" assigns one unique atom  *}
                                          {*per residue  (see <selection>).*}
                 

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



     
                                          {*Use only conformational energy.*}
  flags exclude * include bond angl dihe impr end

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

  for $atom_id in id ( store9 ) loop main
      constraints interaction ( byresidue ( id $atom_id ) ) ( all ) end
      energy end
      vector do ( store1=$ener ) ( id $atom_id )
  end loop main
     
                                   {*Calculate E / rms(E) for each residue.*}
     
  vector show rms (store1) ( store9 )
  vector do (store1=store1/$result) ( store9 )


  {===>} set display=geomplot.print end   {*Write a histogram of energy to *}
                                          {*the file named "geomplot.plot".*}
     
                                                
     
  for $atom_id in id ( store9 ) loop outp
      vector show element (store1) ( id $atom_id )
      if ($result < 1. ) then
         eval ($graph = " " )
      elseif ($result < 2. ) then
           eval ($graph = "X" )
      elseif ($result < 3. ) then
           eval ($graph = "X-X" )
      elseif ($result < 4. ) then
           eval ($graph = "X-X-X" )
      elseif ($result < 5. ) then
           eval ($graph = "X-X-X-X" )
      elseif ($result < 6. ) then
           eval ($graph = "X-X-X-X-X" )
      elseif ($result < 7. ) then
           eval ($graph = "X-X-X-X-X-X" )
      elseif ($result < 8. ) then
           eval ($graph = "X-X-X-X-X-X-X" )
      elseif ($result < 9. ) then
           eval ($graph = "X-X-X-X-X-X-X-X" )
      elseif ($result  < 10. ) then
           eval ($graph = "X-X-X-X-X-X-X-X-X" )
      else
           eval ($graph = "X-X-X-X-X-X-X-X-X-X-->" )
      end if
      vector show elem ( resid ) ( id $atom_id )
      evaluate ( $resid=$result )
      vector show elem ( segid ) ( id $atom_id )
      evaluate ( $segid=$result )
      display  $segid  $resid   $graph
  end loop outp

  {===>} set display=geomplot.list end           {*Write a list file to the*}
                                                 {*specified file.         *} 
  
  evaluate ($number=0)
  for $atom_id in id ( store9 ) loop out2
      vector show element (store1) ( id $atom_id )
      evaluate ( $arg=$result )
      evaluate ($number=$number+1)
      display   $number  $arg 
  end loop out2

  set echo=on end
  set message=on end

  stop