#!/usr/bin/perl
#generate povray input file from DALTON output file
#
# --------------------------------------------------------------------------------
# Revision log:
#
# 20020915 [VB]
# 		Minor modifications, typos, fixed weird aspect ratio
#
# --------------------------------------------------------------------------------
#
#

$file="";
$ffnameprocess="";
$fffile="";
$go=0;
$blocknr=0;
$count=0;
$rotx=0;
$roty=0;
$rotz=0;
$nratoms=0;
$blockcount=0;
$superblock=0;
$onlyfirst=0;
$firstblock=0;
$dist=0;
$size=0;
$makeaxes=0;
@color=0;
$govib=0;
$vibnr=0;
@normmodes=0;
$help1=0;
$help2=0;
$help3=0;
$help4=0;
$vibit=0;
$glass="";
$vibscale="";
$res="";
#____________________________command line__________________________

if (@ARGV[0] eq ""){

    # help message
    print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    print "This perl script ought to generate a molecule.pov POV-Ray \n";
    print "file and a molecule.tga picture from DALTON output files.\n";
    print "Note that the script identifies atoms by their labels \n";
    print "and not by their nuclear charges as DALTON does with \n";
    print "xxxx.mol files. So labels in the xxxx.mol should be given\n";
    print "as:\n";
    print "   H#1       0.00000000    -1.44761450     1.07875057 \n";
    print "   O         0.00000000     0.00000000    -0.13594222 \n";
    print "   H2        0.00000000     1.44761450     1.07875057 \n";
    print "rather than:\n";
    print "   Peter     0.00000000    -1.44761450     1.07875057 \n";
    print "   Paul      0.00000000     0.00000000    -0.13594222 \n";
    print "   Marry     0.00000000     1.44761450     1.07875057 \n";
    print "when you are using dal2pov.\n\n";
    print "Thank you.                                       AAA 1998\n";
    print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    print "\n";
    print "usage: perl dal2pov.perl infile flag1 flag2 \n";
    print "\n";
    print "infile: DALTON input/output file  \n";
    print "flag1 : specifier for infile  \n";
    print "        1 =  geometry optimization output (*OPTIMIZE)\n";
    print "        2 =  final geometry most forms of DALTON output (*WALK etc.) \n";
    print "        3 =  input geometry from DALTON output\n";
    print "        4 =  display path of geometry optimization\n";
    print "        5 =  generate a picture displaying the normal mode\n";
    print "             given as third argument like for normal mode 12:\n";
    print "                perl dal2pov 5 0 12\n";
    print "        6 =  same as 5) but display the vibration in an animation \n";
    print "flag2 : specifier for interactive mode \n";
    print "        1 =  ask for model, animation and axes and size of picture\n";
    print "        2 =  ask for 1), zoom, rotation and scaling of vibrations \n";
    print "        3 =  ask for 1),2) dummy atoms, background color and glass-bonds\n";
    print "             (dummy atoms have to contain \"xx\" in their label !!) \n";
    print "        everything else => use default values\n";
    die("error: no argument\n");
#________________configurating task and error messages__________________

}else{

 $task=@ARGV[1];

    if ($task == ""){
    print "no specifier specified, this probably won't work. . .\n ";
    $task=2;
    };

  if ($task == 1){  $contens="Final geometry"};
  if ($task == 2){  $contens="Molecular geometry"};
  if ($task == 3){  $contens="Total number of coordinates:"};
  if ($task == 4){  $contens="Molecular geometry"};

  if ($task == 6){
     $vibit="y";
     $task = 5;
       };

  if ($task == 5){
     $contens="Molecular geometry";
      if(@ARGV[3]==""){
         print "no mode specified, coosing first one by default . . . \n";
         $vibnr=1;
       }else{
         $vibnr=@ARGV[3];
       };

   };


           if(@ARGV[2]==""){
         print "using default values . . . . \n";
           }else{

#______________________interactive mode for configuration__________________

  if(@ARGV[2]>=1){
print "\nWhat kind of model do you want ?\n";
print "1 = ball and stick (default)\n";
print "2 = sticks\n";
print "3 = space filling\n";
  $mag=<STDIN>;
 chop($mag);

print "animate scene in rotation ? (y/n)\n";
            $moveit=<STDIN>;
            chop($moveit);

print "display axes of output orientatiom ? (y/n)\n";
            $makeaxes=<STDIN>;
            chop($makeaxes);
      };

print "\n";
print "what size shall the Pov. picture have ?\n";
print "1 = 120 \n";
print "2 = 320 \n";
print "3 = 640 \n";
print "4 = 800 \n";
  $res=<STDIN>;
 chop($res);



  if(@ARGV[2]>=2){
    if($moveit=~/y/){
     print "molecule will be animated rotating . . .\n";
    }else{
       print "rotation of molecule x: (-170 for example)";
       $rotx=<STDIN>;
       chop($rotx);
       print "rotation of molecule y: (55 for example)";
       $roty=<STDIN>;
       chop($roty);
       print "rotation of molecule z: (you know...)";
       $rotz=<STDIN>;
       chop($rotz);
    };

print "zoom factor : (10 for H2O, 100 for DNA)\n";
            $zoom=<STDIN>;
            chop($zoom);


  if($task == 5 || $task == 6){
  print "scaling factor for vib. amplitude : (default is o.k. in most cases)\n";
  $vibscale=<STDIN>;
  chop($vibscale);
  };

};

if(@ARGV[2]>=3){
print "display dummy atoms ? (y/n)\n";
            $dum=<STDIN>;
            chop($dum);

print "background color  : \n";
print "type 3 numbers <1 >0 for red/green/blue \n";
print "( \"1.0 1.0 0.0\" for yellow, see colors.inc ) \n";
            $incolor=<STDIN>;
            chop($incolor);

print "make bonds of perspex ( default for 5) ) ? (y/n)\n";
            $glass=<STDIN>;
            chop($glass);
};
print "fishy input, using default . . .";
};

#_______________________open output : molecule.pov _____________________

    # file name
    $file=@ARGV[0];
# print "-----",$task,"\n";
# print "-----",$vibit,"\n";
# print "\n",@ARGV[0],"\n";
# print @ARGV[1],"\n";
# print @ARGV[2],"\n";
# print @ARGV[3],"\n";
    #povray input file
    $fffile = "molecule.pov";
    open(LOG,"$file");

#_______________________reading output__________________

    while (<LOG>){

            #special format for reading .mol file
#           if ($task==&){
#               if ($_ =~ /^[A-Z]....\s*[+-]*[0-9]*.[0-9]\s*/){
#               print "xxxxxxx>",$_;
#               s/\n/ /;
#      	         s/\#/ /;
#                s/\*/ /;
#	        s/\s[0-9]\s//;
#		s/^\s+//;
#                #lines with atom type and xyz go in here
#		$block="$block$_";
#	        };
#
#	     } else {

          #reading .out files

	 if ($go == 1){
            #stop reading if
	    if ($_ =~ /[a-zA-Z][a-zA-Z][a-zA-Z][a-zA-Z][a-zA-Z]/){
          $go=0;
          $onlyfirst=1;
#          print $_;
          };
          };

         #stop reading normal coordinate
	 if ($govib == 1){
            #if
	    if ($_ =~ /[a-zA-Z][a-zA-Z][a-zA-Z][a-zA-Z][a-zA-Z]/){
          $govib=0;
          };
          };

        #reading of normal mode vectors
      	if ($govib == 1){
  # print $_,"\n";
               if ($_ =~ /^...........[x]\s\s/){
                s/\n/ /;
                s/[0-9]*\#[0-9]*/ /;
                s/\*/ /;
      	        s/\s[a-z]\s//;
                s/\s[A-Z]\s//;
	        s/\s[0-9][0-9]*\s//;
                s/\s[0-9]\s//;
     		s/^\s+//;
  #print $_,"\n";
     		$normvibx="$normvibx$_";
	      };

  if ($_ =~ /^...........[y]\s\s/){
                s/\n/ /;
                s/[0-9]*\#[0-9]*/ /;
                s/\*/ /;
      	        s/\s[a-z]\s//;
                s/\s[A-Z]\s//;
	        s/\s[0-9][0-9]*\s//;
                s/\s[0-9]\s//;
     		s/^\s+//;
  #print $_,"\n";
     		$normviby="$normviby$_";
	      };

   if ($_ =~ /^...........[z]\s\s/){
                s/\n/ /;
                s/[0-9]*\#[0-9]*/ /;
                s/\*/ /;
      	        s/\s[a-z]\s//;
                s/\s[A-Z]\s//;
	        s/\s[0-9][0-9]*\s//;
                s/\s[0-9]\s//;
     		s/^\s+//;
  #print $_,"\n";
     		$normvibz="$normvibz$_";
	      };

             };

        #reading of lines with coordinates
	if ($go == 1){

	    if ($task == 3){
              if ($onlyfirst == 0){
              if ($_ =~ /^\s\s\s*[0-9]+\s\s/){
#                print $_,"\n";
                s/\n/ /;
                s/[0-9]*\#[0-9]*/ /;
                s/\*/ /;
      	        s/\s[xyz]\s//;
	        s/\s[0-9][0-9]*\s//;
                s/\s[0-9]\s//;
     		s/^\s+//;
#                print $_,"\n";
                #lines with atom type and xyz go in here
		$firstblock="$firstblock$_";
	      };
	      };
            } else{
#            print $_;
       	      if ($_ =~ / [A-Za-z][a-z]*[a-z]*/){
                s/\n/ /;
	        s/\s[0-9]\s//;
		s/^\s+//;
                #lines with atom type and xyz go in here
		$block="$block$_";
                $superblock="$superblock$_";
	      };
             };
        };

         #look for errors in Input
	 if ($_ =~ /SEVERE ERROR/){
          print"\n WARNING : Error in DALRON Output, this probably won`t work !!!!!\n";
	};

        #start reading if
	 if ($_ =~ /$contens/){
            $go=1;
            $block=0;
	};

        if ($task==5){
        if ($_ =~ /Normal Coordinates/){
#            print $_,"\n";
            $govib=1;
        };
      };
#       };
      };
    close(LOG);

#________________________processing of data____________________


        if ($task == 3){
	@blocklines=split(/\s+/,$firstblock);
        }else{
        #split in-srting into array
	@blocklines=split(/\s+/,$block);
	@blocklines2=split(/\s+/,$superblock);
        $numbersteps=($#blocklines2+1)/($#blocklines+1);
        };

	@color=split(/\s+/,$incolor);

       #determine the farest point from 0,0,0
         $k=1;
         for ($j=1;$j<(($#blocklines+1)*0.25)+1;$j++){
         $dist=sqrt((@blocklines[$k]**2)+(@blocklines[$k+1]**2)+(@blocklines[$k+2]**2));
         if($dist>$size){$size=$dist};
         $k=$k+4;
         };


         if ($task == 5){

	@normmodesx=split(/\s+/,$normvibx);
	@normmodesy=split(/\s+/,$normviby);
	@normmodesz=split(/\s+/,$normvibz);


        $help1=$vibnr;
	  while($help1>0){

        $help1=$help1-5;
        $help2=(1+$help2);
        };

        $help3=(($#blocklines+1)*0.75)-6;
	  while($help3>0){
         $help3=$help3-5;
        $help4=1+$help4;
        };

          $help5=5;
	  if ($help4 eq $help2){
           $help5=(5+$help3)};

#        print$help1,"..",$help2,"\n";
#        print$help3,"..",$help4,"\n";
#        print"..",$help5,"\n";
        $ampli=0;

        $k=0;
        $m=0;
        for ($j=1;$j<(($#blocklines+1)*0.25)+1;$j++){
        @normvec[$k]  =@normmodesx[$m+4+$help1+(($#blocklines+1)*0.25*$help5*($help2-1))];
        @normvec[$k+1]=@normmodesy[$m+4+$help1+(($#blocklines+1)*0.25*$help5*($help2-1))];
        @normvec[$k+2]=@normmodesz[$m+4+$help1+(($#blocklines+1)*0.25*$help5*($help2-1))];


         if($ampli < (sqrt(@normvec[$k]**2+@normvec[$k+1]**2+@normvec[$k+2]**2))){
            $ampli = (sqrt(@normvec[$k]**2+@normvec[$k+1]**2+@normvec[$k+2]**2));
	  };

        $k=$k+3;
        $m=$help5+$m;
#        print"....",$k," ",$j;
         };
         };

 if ($vibit=~/y/){$task=6};

#print "___",$#blocklines2+1,"  ",$#blocklines+1,"   ",$numbersteps,"\n";
#print ;
#print  " --->",$#blocklines,"\n";
#for ($j=0;$j<30;$j++){
#print"-----\n";
#print@normvec[$j];
#print "\n";
#};

#____________________________defining default values____________

#background color
if (@color[0]==""){@color[0]=0.2};
if (@color[1]==""){@color[1]=0.2};
if (@color[2]==""){@color[2]=0.4};

#rotation
if ($rotx==""){$rotx=0};
if ($roty==""){$roty=0};
if ($rotz==""){$rotz=0};

#space fill, stick and/or ball
if ($mag==2){
$mag=0.2;
	   }else{
if ($mag==3){
$mag=3;
	   }else{
$mag = 1.15;
};
};

#zoom factor
if ($zoom==""){$zoom=($size*3.5)};

#scale factor for normalmodes
if (! $vibscale == ""){$ampli = $ampli/$vibscale};

#dummy atoms
if ($dum==""){
    $dum=0;
    };

#glass bonds
if ($task == "5" && $glass == ""){
     $glass="y";
    }else{;
    $glass=0;
};


#_______________________output generation_______________________

 open(FF,">$fffile");

         #povray header
 print  FF " // Persistence Of Vision raytracer version 3.0 sample file.   \n";
 print FF " #version 3.0                                        \n";
 print FF " \n";
 print  FF "global_settings { assumed_gamma 2.2 }             \n";
 print FF " \n";
 print  FF "#include \"colors.inc\"     // file for colors \n";
 print FF  "#include \"shapes.inc\"     // file for sphere, cylinder etc.   \n";
 print  FF "#include \"textures.inc\"   // take a look at the file textures.v2  \n";
 print  FF "#include \"chars.inc\"      // file for letters and numbers     \n";
 print FF " \n";
 print FF "// this is the camera part, to change the zoom edit the z value of location \n";
 print FF "camera {                        \n";
 print FF  "    location  <0, 0,-",$zoom,">\n";
 print FF "    direction <0, 0,  1>        \n";
 print FF "    up        <0,  1,  0>       \n";
 print FF "    right     x*1.3333          \n";
 print FF "    look_at   <0, 0, 0>         \n";
 print FF "    }                           \n";
 print FF " \n";
 print FF "background {color red ",@color[0]," green ",@color[1]," blue ",@color[2]," }    \n";
  print FF " \n";
print  FF "light_source {<20, 30, -100> colour 1}         \n";
 print FF " \n";
 print FF " \n";
if ($moveit=~/y/){

 print FF "// these are the clock values for a rotation animation \n";
 print  FF "#declare Clock360 = 360*clock    // normalize steps \n";
 print  FF "#declare ClockRotz = Clock360*z   \n";
 print  FF "#declare ClockRotx = Clock360*x   \n";
 print  FF "#declare ClockRoty = Clock360*y   \n";
  print FF " \n";
                };
#write variables with coordinates time dependent on vibration
if ($vibit=~/y/){
 print FF "// these are the clock values for an animated vibration \n";
 print  FF "#declare ClockPI = 2*3.1415926*clock  // clock runs over 2 Pi \n";
 print  FF "#declare ClockVib = sin(ClockPI)      // sine amplitude \n";
 print FF " \n";
  $k=0;
  $l=0;

 print FF "// these are the coordinates of the atoms plus scaled cartesian vibrational eigenvectors of \n";
 print FF "// multiplied by the sine amplitude from the clockset \n";
  for ($j=1;$j<(($#blocklines+1)*0.25)+1;$j++){
#print" .   ",@normvec[$k-$l],"  .  ",@normvec[$k+1-$l],"  .  ",@normvec[$k+2-$l],"\n";
 print FF " \n";
   print FF "#declare x",$j,"=",@blocklines[$k+1],"+",(0.5/$ampli)*@normvec[$k-$l],"\*ClockVib","\n";
   print FF "#declare y",$j,"=",@blocklines[$k+2],"+",(0.5/$ampli)*@normvec[$k+1-$l],"\*ClockVib","\n";
   print FF "#declare z",$j,"=",@blocklines[$k+3],"+",(0.5/$ampli)*@normvec[$k+2-$l],"\*ClockVib","\n";
 print FF " \n";
    $l=$l+1;
    $k=$k+4;
   };
}else{

#write static variables with coordinates
 print FF "// cartesian coordinates of the atoms \n";
$k=1;
for ($j=1;$j<(($#blocklines+1)*0.25)+1;$j++){
 print FF " \n";
print FF "#declare x",$j,"=",@blocklines[$k],"\n";
print FF "#declare y",$j,"=",@blocklines[$k+1],"\n";
print FF "#declare z",$j,"=",@blocklines[$k+2],"\n";

$k=$k+4;
};
};

 print FF " \n";
 print FF " \n";
 print  FF "union{ // this part is treated as a union \n";

#______display path of geometry optimization with yellow lines_______________________

 if ($task == 4){

$k=1;
for ($j=1;$j<((($#blocklines+1)*0.25)*($numbersteps-1));$j++){
#print $j,"   :  ";

  #treshold to avoid degenerate cylinders
  if(sqrt((@blocklines2[$k]-@blocklines2[$k+($#blocklines+1)])**2+(@blocklines2[$k+1]-@blocklines2[$k+1+($#blocklines+1)])**2+(@blocklines2[$k+2]-@blocklines2[$k+2+($#blocklines+1)])**2)>0.00001){
 print FF " \n";
print FF "cylinder{<",@blocklines2[$k],",",@blocklines2[$k+1],",",@blocklines2[$k+2]," >, \n";
print FF  "<",@blocklines2[$k+($#blocklines+1)],",",@blocklines2[$k+1+($#blocklines+1)],",",@blocklines2[$k+2+($#blocklines+1)]," >,0.03 \n";
print FF " pigment { color red 1 green 1 blue 0 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
};
$k=$k+4;
};
};

 print FF " \n";
 print FF " \n";
 print FF "// these are the spheres for the atoms working like this : \n";
 print FF "// sphere     x,y,z  and radius \n";
 print FF "// pigment    color of the sphere  \n";
 print FF "// finish     properties of the sphere's surface\n";
  $k=0;
  for ($j=1;$j<(($#blocklines+1)*0.25)+1;$j++){

#dummy atoms, check it out
if(@blocklines[$k]=~/xx/){
if($dum=~/y/){
  print FF " \n";
print FF "cylinder { <x",$j,"-0.01,y",$j,"-0.01,z",$j,"-0.01>, <x",$j,"+0.01,y",$j,"+0.01,z",$j,"+0.01>,", 0.3,"\n";
print FF " pigment { color red 1.0 green 0.6 blue 0.2 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
 print FF " \n";
print FF "cylinder { <x",$j,"-0.01,y",$j,"-0.01,z",$j,"-0.01>, <x",$j,"+0.01,y",$j,"+0.01,z",$j,"+0.01>,",0.3,"\n";
print FF " pigment { color red 1.0 green 0.6 blue 0.2 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }\n";
print FF "rotate <90,0,0>}","\n";
 print FF " \n";
print FF "cylinder { <x",$j,"-0.01,y",$j,"-0.01,z",$j,"-0.01>, <x",$j,"+0.01,y",$j,"+0.01,z",$j,"+0.01>,",0.3,"\n";
print FF " pigment { color red 1.0 green 0.6 blue 0.2 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }\n";
print FF "rotate <0,0,90>}","\n";
 print FF " \n";
};

}else{

#other atoms
    if(@blocklines[$k]=~/H/){
 print FF " \n";
print FF "sphere { <x",$j,",y",$j,",z",$j," >,",0.4*$mag,"\n";
print FF " pigment { color red 1 green 1 blue 1 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
    }else{

if(@blocklines[$k]=~/O/){
 print FF " \n";
print FF "sphere { <x",$j,",y",$j,",z",$j," >,",0.6*$mag,"\n";
print FF " pigment { color red 1 green 0 blue 0 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
}else{

if(@blocklines[$k]=~/N/){
 print FF " \n";
print FF "sphere { <x",$j,",y",$j,",z",$j," >,",0.6*$mag,"\n";
print FF " pigment { color red 0 green 0 blue 1 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
}else{


if(@blocklines[$k]=~/C/){

  if(@blocklines[$k]=~/Cl/){
 print FF " \n";
  print FF "sphere { <x",$j,",y",$j,",z",$j," >,",0.8*$mag,"\n";
  print FF " pigment { color red 0 green 1 blue 0 }","\n";
  print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
  }else{
 print FF " \n";
print FF "sphere { <x",$j,",y",$j,",z",$j," >, ",0.7*$mag,"\n";
print FF " pigment { color red 0.1 green 0.1 blue 0.1 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
  };

}else{

if(@blocklines[$k]=~/Br/){
 print FF " \n";
print FF "sphere { <x",$j,",y",$j,",z",$j," >,",1.0*$mag,"\n";
print FF " pigment { color red 0.4 green 0.4 blue 0.6 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
}else{

if(@blocklines[$k]=~/F/){
 print FF " \n";
print FF "sphere { <x",$j,",y",$j,",z",$j," >,",0.5*$mag,"\n";
print FF " pigment { color red 0.2 green 1 blue 0 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
}else{



if(@blocklines[$k]=~/S/){
 print FF " \n";
print FF "sphere { <x",$j,",y",$j,",z",$j," >,",0.8*$mag,"\n";
print FF " pigment { color red 1 green 1 blue 0 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";

}else{
#undefined, just grey
 print FF " \n";
print FF "sphere { <x",$j,",y",$j,",z",$j," >,",0.5*$mag,"\n";
print FF " pigment { color red 0.5 green 0.5 blue 0.5 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
};};};};};};};};

$k=$k+4;
};

 print FF " \n";
 print FF " \n";
 print FF "// these are the cylinders for the bonds working like this : \n";
 print FF "// cylinder   x,y,z of cylinder startpoint, x,y,z of end and radius \n";
 print FF "// finish     properties of the sphere's surface\n";
 print FF "// pigment    color  \n";
 print FF "// in case of glass bonds the finish and ambient are replaced by  \n";
 print FF "// the texture statement followed by the texture identifier Glass3 \n";
#now for the bonds, everything closer than 2.8 au
$k=1;
for ($j=1;$j<(($#blocklines+1)*0.25)+1;$j++){
$m=$k+4;
for ($l=$j+1;$l<(($#blocklines+1)*0.25)+1;$l++){

$dist=sqrt((@blocklines[$k]-@blocklines[$m])**2+(@blocklines[$k+1]-@blocklines[$m+1])**2+(@blocklines[$k+2]-@blocklines[$m+2])**2);

if ($dist<2.8){
 print FF " \n";
print FF "cylinder { <x",$j,",y",$j,",z",$j,">, <x",$l,",y",$l,",z",$l,">, 0.15","\n";
if ($glass=~/y/){
print FF " texture {Glass3}","\n";
}else{
print FF " finish {","\n";
print FF " ambient 0.1 ","\n";
print FF " diffuse 0.6","\n";
print FF "        } ","\n";
print FF " pigment {color red 0.99 green 0.99 blue 1.0} ","\n";
};
print FF "     }","\n";
};

$m=$m+4;
};
$k=$k+4;
};


#here come the axes
 if($makeaxes =~/y/){
  print FF " \n";
print FF " \n";
print FF "// this is the coordinate system of the DALTON Geometry \n";
print FF "// the x,y and z axis are displayed as follows :\n";
print  FF "union{ \n";
print FF " \n";
print FF "#declare e =",$size+0.1*$size,"// this is somehow scale of the axis-size \n";
print FF " \n";
 print FF "// each axis is a cylinder with a cone on top : \n";
 print FF "// X - axis--------------------------- \n";
print FF "cone {  \n";
print FF "  <e,0,0>, e/36 // cone origin and basewidth \n";
print FF "  <e+(e/9),0,0>, 0.0  // cone tip and tipwidth \n";
print FF " pigment {Silver}\n";
print FF "}\n";
print FF " \n";
  print FF "// cylinders just the same as for the bonds \n";
print FF " cylinder { <0,0,0>, <e,0,0>, e/60\n";
print FF "  finish {\n";
print FF "  ambient 0.1 \n";
print FF "  diffuse 0.6\n";
print FF "  } \n";
print FF "  pigment {Silver} \n";
print FF "   }\n";
print FF " \n";

 print FF "// Y - axis--------------------------- \n";
print FF "cone {  \n";
print FF "  <0,e,0>, e/36\n";
print FF "  <0,e+(e/9),0>, 0.0\n";
print FF " pigment {Silver}\n";
print FF "}\n";
print FF " \n";

print FF " cylinder { <0,0,0>, <0,e,0>, e/60\n";
print FF "  finish {\n";
print FF "  ambient 0.1 \n";
print FF "  diffuse 0.6\n";
print FF "  } \n";
print FF "  pigment {Silver} \n";
print FF "   }\n";
print FF " \n";

print FF "// Z - axis--------------------------- \n";
print FF "cone {  \n";
print FF "  <0,0,e>, e/36\n";
print FF "  <0,0,e+(e/9)>, 0.0\n";
print FF " pigment {Silver}\n";
print FF "}\n";
print FF " \n";

print FF " cylinder { <0,0,0>, <0,0,e>, e/60\n";
print FF "  finish {\n";
print FF "  ambient 0.1 \n";
print FF "  diffuse 0.6\n";
print FF "  } \n";
print FF "  pigment {Silver} \n";
print FF "   }\n";
print FF " \n";

print FF "// a sphere at the origin to connect the cylinders \n";
print FF "sphere { <0,0,0>, e/60;";
print FF "  pigment {Silver} \n";
print FF " finish { ambient 0.1 diffuse 0.6  }}","\n";
print FF "   }\n";
print FF " \n";
print FF " \n";

print FF "// these are the characters \n";
print FF " \n";
print FF "union{ \n";
print FF "          char_Z \n";
print FF "     translate <0,-2.5,-0.5> \n";
print FF "          pigment{ color Yellow } \n";
print FF "scale <e*0.02,e*0.02,e*0.02> // they have to be scaled \n";
print FF "translate <0,0,e+(e/6)>      // and translated above the tip of the cones  \n";
print FF "        } \n";
print FF " \n";

print FF "union{ \n";
print FF "          char_Y \n";
print FF "     translate <0,-2.5,-0.5> \n";
print FF "          pigment{ color Yellow } \n";
print FF "scale <e*0.02,e*0.02,e*0.02> \n";
print FF "translate <0,e+(e/6),0> \n";
print FF "        } \n";
print FF " \n";

print FF "union{ \n";
print FF "          char_X \n";
print FF "     translate <0,-2.5,-0.5> \n";
print FF "          pigment{ color Yellow } \n";
print FF "scale <e*0.02,e*0.02,e*0.02> \n";
print FF "translate <e+(e/6),0,0> \n";
print FF "        } \n";
print FF " \n";
print FF "//---------------------- \n";

};
 #vectors to display normal modes
 if ($task==5){
  $k=0;
  $l=0;
  for ($j=1;$j<(($#blocklines+1)*0.25)+1;$j++){

  # treshold to avoid degenerate cylinders

  if(sqrt(((@normvec[$k-$l])**2)+((@normvec[$k+1-$l])**2)+((@normvec[$k+2-$l])**2))>0.00001){
print FF " \n";
print FF "cylinder{<",@blocklines[$k+1],",",@blocklines[$k+2],",",@blocklines[$k+3]," >, \n";
print FF "<",@blocklines[$k+1]+@normvec[$k-$l]*(1/$ampli),",",@blocklines[$k+2]+@normvec[$k+1-$l]*(1/$ampli),",",@blocklines[$k+3]+@normvec[$k+2-$l]*(1/$ampli)," >,0.03 \n";
print FF "pigment { color red 1 green 1 blue 0 }","\n";
print FF " finish { ambient 0.3 diffuse 0.7 phong 1 }}","\n";
print FF " \n";
print FF "cone {<",@blocklines[$k+1]+@normvec[$k-$l]*(1/$ampli),",",@blocklines[$k+2]+@normvec[$k+1-$l]*(1/$ampli),",",@blocklines[$k+3]+@normvec[$k+2-$l]*(1/$ampli)," >,0.05 \n";
print FF "<",(@blocklines[$k+1]+@normvec[$k-$l]*(1/$ampli)*1.2),",",(@blocklines[$k+2]+@normvec[$k+1-$l]*(1/$ampli)*1.2),",",(@blocklines[$k+3]+@normvec[$k+2-$l]*(1/$ampli)*1.2)," >,0.0 \n";
print FF " pigment {Yellow}\n";
print FF "}\n";
print FF " \n";
};
$k=$k+4;
$l=$l+1;
  };
 };

if ($moveit=~/y/){
#print  FF " rotate ClockRotz  \n";
print FF " \n";
print FF " \n";
print  FF " rotate ClockRotx \n";
print  FF " rotate 2*ClockRoty \n";
   }else{
print FF " \n";
print  FF " rotate <",   $rotx,    ",",     $roty,   ",",     $rotz,    "> \n";
   };

print FF "   }\n";
print FF " \n";

   close(FF);


#__________________call povray : molecule.pov -> molecule.tga_________

 if ($vibit=~/y/){$moveit="y"};

if ($moveit=~/y/){
   if ($res eq ""){
    system('povray +Q9 rotate.ini res320.ini -I molecule.pov');
   }else{
        if ($res=~/1/){    system('povray +Q9 rotate.ini res120.ini -I molecule.pov')};
        if ($res=~/2/){    system('povray +Q9 rotate.ini res320.ini -I molecule.pov')};
        if ($res=~/3/){    system('povray +Q9 rotate.ini res640.ini -I molecule.pov')};
        if ($res=~/4/){    system('povray +Q9 rotate.ini res800.ini -I molecule.pov')};
        };
}else{
   if ($res eq ""){
    system('povray res640.ini -I molecule.pov');
   }else{
       if ($res=~/1/){  system('povray res120.ini -I molecule.pov')};
       if ($res=~/2/){ system('povray res320.ini -I molecule.pov')};
       if ($res=~/3/){ system('povray res640.ini -I molecule.pov')};
       if ($res=~/4/){ system('povray res800.ini -I molecule.pov')};
       };
};
die("dal2pov finished\n");
}

#declare s1 = 0
#declare s2 = 0
#declare s3 = 0
#declare e1 = 1
#declare e2 = 1.3
#declare e3 = 0.6

#cone {
#  <e1,e2,e3>, 0.2
#  <e1+((e1-s1)/4),e2+((e2-s2)/4),e3+((e3-s3)/4)>, 0.0
# pigment {NeonBlue}
#   translate <0,0,0>
#}
# cylinder { <s1,s2,s3>, <e1,e2,e3>, 0.1
#  finish {
#  ambient 0.1
#  diffuse 0.6
#  }
#  pigment {NeonBlue}
#   }