CCL:G: Optimization with imaginary frequency

[René Kanters <rkanters=richmond.edu>]

 Sent to CCL by: =?ISO-8859-1?Q?René_Kanters?= [rkanters++richmond.edu]
 Hi Flavio,
I hacked together a little perl scrip that will look at the first set
 frequency information in a gaussian calculation and generate a new
 set of xyz coordinates in which the structure is moved along a
 vibrational coordinate for a certain amount (see the command line
 arguments in the script below).
Note that it may be necessary, depending on the symmetry of the
 saddlepoint you got to create a new input with going in the positive
 and another one in the negative direction of the imaginary frequency
 since you can't be sure which one will lead you to a lower minimum.
 I hope this helps.
 René
PS it's not pretty (I am not a perl programmer), but it has worked
 for me.
 #! /usr/bin/perl
 ###############################################
 # script to take a gaussian log file output and generate structure
 # information (cartesian) for a new structure with a displacement
 # along a frequency coordinate in the direction along the '-v' vibration
 # for a displacement of '-d' times the vibrational vector
 #
 # command line arguments:
 #
 # -h        : prints help and quits.
 # -d float  : displacement along vibrational mode default 0.1
 # -v int    : vibration number to use default 1
 # required : filename (only sees last one)
 ################################################
 #
 $test = 0;  # for debugging purposes (prints out extra stuff)
 # set up the defaults
 $vibration = 1;
 $displacement = 0.1;
 $file = "";
 # 'constants' not using the 'use constant NAME => value;' method.
 $C_STRUCTURE = "Standard orientation:";
 $C_ENDSTRUCTURE = "---------";
 $C_FREQUENCY = "Frequencies --";
 while ($#ARGV >= 0) {
   if ($ARGV[0] eq "-d")    { $displacement = $ARGV[1]; shift |-|ARGV;
 }
   elsif ($ARGV[0] eq "-v") { $vibration = $ARGV[1]; shift |-|ARGV; }
   elsif ($ARGV[0] eq "-h") { &usage(); exit;}
   else { $file = $ARGV[0];}
   shift |-|ARGV;
 }
 if ($file eq "") { &usage(); exit;}
 open(IN,$file) || die "Could not open file '$file'.\n";
 while (1) {
   # read until I get a structure (should happen first) or a frequency
  until ($_ = <IN>, $_ =~ /$C_STRUCTURE/ or $_ =~ /$C_FREQUENCY/
 or
 eof) { }
   if ($_ =~ /$C_STRUCTURE/) {  # READ THE STRUCTURE
    $nat = -1;   # reset the atom index counter : only need last
 structure
     $at = ();    # clear out the arrays of atomnumber, x, y and z
     $x = ();
     $y = ();
     $z = ();
     skiplines(4); # position read line to first atom in the list
     while ($_ = <IN> and not $_ =~ /$C_ENDSTRUCTURE/ ) {
       |-|tokens = split(' ', $_);
       $nat++;
       $at[$nat] = $tokens[1];
      $xind = $#tokens - 2; # needed if the atomic type is missing
 (happens sometimes).
       $x[$nat] = $tokens[$xind];
       $y[$nat] = $tokens[$xind+1];
       $z[$nat] = $tokens[$xind+2];
     }
   } elsif ($_ =~ /$C_FREQUENCY/) { # READ THE FREQUENCY INFORMATION
     # since there are three vibrations per line, I can determine how
    # many of these I need to skip before I can read the frequency I
 need.
     $skipsize = 4 + ($nat+1) + 3;  # skip rest header, vibrations
 and start header
     for ($i=0; $i<int($vibration/3); $i++) { skiplines($skipsize); }
     skiplines(4);
    $xind = 2 + (($vibration-1) % 3) * 3; # position in the tokens
 of vib vector
     if ($test eq 1) { print "vibration = $vibration, xind = $xind\n";
 }
     for ($i=0;$i<=$nat;$i++) {
       $_ = <IN>;
       |-|tokens = split(' ',$_);
       if ($test eq 1 ) {
         print "input $_";
        print "($x[$i],$y[$i],$z[$i]) + $displacement x ($tokens
 [$xind],$tokens[$xind+1],$tokens[$xind+2])\n";
       }
      printf "%4d  % .6f  % .6f  % .6f\n",$at[$i], $x[$i]+
 $displacement*$tokens[$xind], $y[$i]+$displacement*$tokens[$xind+1],
 $z[$i]+$displacement*$tokens[$xind+2];
     }
     exit; # AND WE ARE DONE!!!!
   } elsif (eof) { close IN; exit; }
 }
 ##################################################
 # skiplines : argument n: number of lines to skip
 ##################################################
 sub skiplines {
   local $n, $i;
   $n = shift(|-|_);
   for ($i=0;$i<$n;$i++) { <IN> || die "skiplines encountered
 EOF."; }
 }
 ##################################################
 # usage prints a message about how to use sp2g03
 ##################################################
 sub usage {
 print <<ENDUSAGE;
 Usage: nfreq [-h] [-d displacement] [-v vibrationnnumber] filename
 -d : displacement step along vibrational mode default 0.1
 -v : index of vibrational mode to follow default 1
 -h : prints this help
Will create the structure '-d' along the '-v'th vibrational mode of
 the stucture.
 ENDUSAGE
 }
On Jul 23, 2008, at 3:39 PM, flavio forti fforti.:.fortisrl.com.ar
 wrote:
>  Sent to CCL by: "flavio  forti" [fforti^_^fortisrl.com.ar]
> Hi everybody. I am trying to optimize a large set of histamine
>  monocation  conformers in solution (PCM) with B3LYP/6-31G*. I had
>  some which converged to a saddle point (1 imaginary frequency) but
>  finally, after several reoptimizations with CalcFC they arrived to
>  a minimum. But there is one still resisting. CalcAll would be to
>  expensive, what do you suggest?
>  thanks
>  flavio
> -= This is automatically added to each message by the mailing
>  script =-
>  To recover the email address of the author of the message, please
>  change> Conferences: http://server.ccl.net/chemistry/announcements/
>  conferences/
> Search Messages: http://www.ccl.net/htdig  (login: ccl, Password:
>  search)>