From owner-chemistry@ccl.net Thu Jul 24 13:16:00 2008
From: "=?ISO-8859-1?Q?Ren=E9_Kanters?= rkanters(_)richmond.edu" <owner-chemistry#server.ccl.net>
To: CCL
Subject: CCL:G: Optimization with imaginary frequency
Message-Id: <-37410-080724123206-1131-Yz55mH/NEoc1dA2geEAyXg#server.ccl.net>
X-Original-From: =?ISO-8859-1?Q?Ren=E9_Kanters?= <rkanters=richmond.edu>
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain; charset=ISO-8859-1; delsp=yes; format=flowed
Date: Thu, 24 Jul 2008 12:02:17 -0400
Mime-Version: 1.0 (Apple Message framework v753.1)


Sent to CCL by: =?ISO-8859-1?Q?Ren=E9_Kanters?= [rkanters++richmond.edu]
Hi Flavio,

I hacked together a little perl scrip that will look at the first set =20=

frequency information in a gaussian calculation and generate a new =20
set of xyz coordinates in which the structure is moved along a =20
vibrational coordinate for a certain amount (see the command line =20
arguments in the script below).

Note that it may be necessary, depending on the symmetry of the =20
saddlepoint you got to create a new input with going in the positive =20
and another one in the negative direction of the imaginary frequency =20
since you can't be sure which one will lead you to a lower minimum.

I hope this helps.

Ren=E9

PS it's not pretty (I am not a perl programmer), but it has worked =20
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 =3D 0;  # for debugging purposes (prints out extra stuff)

# set up the defaults
$vibration =3D 1;
$displacement =3D 0.1;
$file =3D "";

# 'constants' not using the 'use constant NAME =3D> value;' method.
$C_STRUCTURE =3D "Standard orientation:";
$C_ENDSTRUCTURE =3D "---------";
$C_FREQUENCY =3D "Frequencies --";

while ($#ARGV >=3D 0) {
   if ($ARGV[0] eq "-d")    { $displacement =3D $ARGV[1]; shift |-|ARGV; }
   elsif ($ARGV[0] eq "-v") { $vibration =3D $ARGV[1]; shift |-|ARGV; }
   elsif ($ARGV[0] eq "-h") { &usage(); exit;}
   else { $file =3D $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 ($_ =3D <IN>, $_ =3D~ /$C_STRUCTURE/ or $_ =3D~ /$C_FREQUENCY/ =
or =20
eof) { }
   if ($_ =3D~ /$C_STRUCTURE/) {  # READ THE STRUCTURE
     $nat =3D -1;   # reset the atom index counter : only need last =20
structure
     $at =3D ();    # clear out the arrays of atomnumber, x, y and z
     $x =3D ();
     $y =3D ();
     $z =3D ();
     skiplines(4); # position read line to first atom in the list
     while ($_ =3D <IN> and not $_ =3D~ /$C_ENDSTRUCTURE/ ) {
       |-|tokens =3D split(' ', $_);
       $nat++;
       $at[$nat] =3D $tokens[1];
       $xind =3D $#tokens - 2; # needed if the atomic type is missing =20=

(happens sometimes).
       $x[$nat] =3D $tokens[$xind];
       $y[$nat] =3D $tokens[$xind+1];
       $z[$nat] =3D $tokens[$xind+2];
     }
   } elsif ($_ =3D~ /$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 =20=

need.
     $skipsize =3D 4 + ($nat+1) + 3;  # skip rest header, vibrations =20
and start header
     for ($i=3D0; $i<int($vibration/3); $i++) { skiplines($skipsize); }

     skiplines(4);
     $xind =3D 2 + (($vibration-1) % 3) * 3; # position in the tokens =20=

of vib vector
     if ($test eq 1) { print "vibration =3D $vibration, xind =3D =
$xind\n"; }
     for ($i=3D0;$i<=3D$nat;$i++) {
       $_ =3D <IN>;
       |-|tokens =3D split(' ',$_);
       if ($test eq 1 ) {
         print "input $_";
         print "($x[$i],$y[$i],$z[$i]) + $displacement x ($tokens=20
[$xind],$tokens[$xind+1],$tokens[$xind+2])\n";
       }
       printf "%4d  % .6f  % .6f  % .6f\n",$at[$i], $x[$i]+=20
$displacement*$tokens[$xind], $y[$i]+$displacement*$tokens[$xind+1], =20
$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 =3D shift(|-|_);
   for ($i=3D0;$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 =20
the stucture.
ENDUSAGE
}

On Jul 23, 2008, at 3:39 PM, flavio forti fforti.:.fortisrl.com.ar =20
wrote:

>
> Sent to CCL by: "flavio  forti" [fforti^_^fortisrl.com.ar]
> Hi everybody. I am trying to optimize a large set of histamine =20
> monocation  conformers in solution (PCM) with B3LYP/6-31G*. I had =20
> some which converged to a saddle point (1 imaginary frequency) but =20
> finally, after several reoptimizations with CalcFC they arrived to =20
> a minimum. But there is one still resisting. CalcAll would be to =20
> expensive, what do you suggest?
> thanks
> flavio
>
>
>
> -=3D This is automatically added to each message by the mailing =20
> script =3D-
> To recover the email address of the author of the message, please =20
> change> Conferences: http://server.ccl.net/chemistry/announcements/=20
> conferences/
>
> Search Messages: http://www.ccl.net/htdig  (login: ccl, Password: =20
> search)>
>