TURBOMOLE Users Forum

General FAQ => General FAQ => Topic started by: uwe on April 30, 2008, 04:29:23 PM

Title: Constrained optimization / scan along internal coordinates
Post by: uwe on April 30, 2008, 04:29:23 PM

Constrained optimization / scan along internal coordinates

To fix internal coordinates and change the value of fixed internal coordinates to do a potential energy surfaces (PES)  scan in TURBOMOLE, one can use either Pulay's internal coordinates (iaut in define), or the redundant internal coordinates (ired in define).

To fix an internal coordinate:

• call define
• read in your coordinates
• Now there are two possibilities to fix internal coordinates:
  
  • For both kinds of internal coordinates (Pulay and redundant):
    
    • enter
         idef
      in the geometry menu of define.
    • define the internal coordinate you want to be fixed by hand (see tutorial for details), .e.g.
         f stre 1 2
      f stands for fixed, stre for strech (bend for bond angles, tors for torsions, outp for out of plane angles, etc.).
      define prints on the screen the names of the internal coordinates:
      
      ENTER INTERNAL COORDINATE DEFINITION COMMAND
      
                       <x> <type> <indices>
      
       WHERE    <x>       = k    f    d    i
                <type>    = stre invr bend outp tors linc linp
                            comp ring pyrm bipy pris cube octa
      
       THESE COMMANDS WILL BE EXPLAINED IN DETAIL IF YOU ENTER  <x> <type> ?  FOR
       SOME CHOICE OF  <x>  AND  <type> , E.G.  k stre ?
      
      As described by define, to get a help about the numbering scheme for the atoms, enter e.g.
        f outp ?
       
    • Go back to the geometry menu (hit a few times <Enter>), and say either:
        ired
      generates internal redundant coordinates taking the fixed internal ones and let define automatically generate a full set of internal redundant coordinates, or
        iaut
       imet
      to generate Pulay's internal coordinates with iaut and check if they are linearly independent with imet.
  • For non-redundant internal coordinates:
    
    • go to the i and then iaut submenu and accept the default settings with *
    • enter imet to check if the number of linearly independent coordinates equals the number of degrees of freedom:
      Look at the first line of the i submenu:
      
      INTERNAL COORDINATE MENU  ( #ideg=45      #k=43   #f=2    #d=2    #i=0     )
      
       imet <a>  : PROVIDE B-MATRIX FOR ACTIVE INTERNAL COORDINATES
                   (CHECK COMPLETENESS AND NUMERICAL QUALITY
                   AND CHANGE REDUNDANT INTERNALS TO display)
       idef      : SUB-MENU FOR INTERACTIVE DEFINITION OF INTERNAL COORDINATES
       ideg <a>  : OUTPUT NUMBER OF TOT. SYMMETRIC INTERNAL DEGREES OF FREEDOM
      [...]
      
      
      #ideg is the number of degrees of freedom, #k is the number of non-fixed internal coordinates and #f the number of fixed internal coordinates. #d and #i are not used.
      
      Hence, you have to make sure that after the imet step #idef = #k + #f
      
      
    • look at the list of internal coordinates with disi, change the status of the internal coordinates with ic, or:
    • run to the end of define, edit the coord file and search for the requested internal coordinate,
        change it and fix it by replacing the k type to an f type, e.g. water might look like this:
         1 f  1.0000000000000 stre    1    2           val=   1.6
         2 k  1.0000000000000 stre    1    3           val=   1.85256
         3 k  1.0000000000000 bend    3    2    1      val= 104.18231

Scan along one internal coordinate:

• first define one or several fixed internal coordinates as described above
• generate several subdirectories, each one representing a geometry at one certain fixed value of the internal coordinate you would like to scan along
• copy your input files to all directories
• in each directory,  change the value of the fixed internal coordinate in the coord file (example for internal redundant coordinates):
  
  $intdef
  # definitions of internal coordinates
     1 f  1.0000000000000 stre    1    2           val=   2.91525
  
  The value of bonds is given in atomic units, those of angles and torsions in degrees. So simply change the value of the internal coordinate, e.g.:
  
     1 f  1.0000000000000 stre    1    2           val=   3.1
  
  If Pulay coordinates are used, the $intdef section of the coord file will also contain all other internal non-fixed coordinates, while in case of internal redundant coordinates the complete set is written to $redundant. Please do not change the values in $redundant, just in $intdef !
  
• run define in each of the subdirectories
• There are again two possibilities now:
  
  
  • in case of Pulay internal coordinates, define will recognize that cartesian and internal coordinates do not match any more, and you will be asked what to use:
    
       CARTESIAN COORDINATES AND VALUES OF INTERNAL COORDINATES DO  N O T   AGREE !
       ENTER COMMAND :
       c  : USE CARTESIAN COORDINATES TO ESTABLISH GEOMETRY (DEFAULT)
       i  : USE VALUES OF INTERNAL COORDINATES TO RE-SHAPE GEOMETRY
    
    If you say i here, the cartesian coordinates will be constructed from the set of internal ones, taking care of the new values of the fixed internal coordinates.
    
    Please note: If you do a geometry optimization with internal coordinates, the values of the internal coordinates in the coord file are not changed! TURBOMOLE will only use the definition of the internal coordinates, not their values!
    So saying i in this part of define will always give you the cartesian coordinates of the structure you had when you have defined the internal coordinates with iaut -> this is different from what you have in cartesian coordinates in the coord file at that point.
    
    
  • In case of internal redundant coordinates:
    
       call ired in the geometry menu which will read in the fixed coordinates and build the remaining internal redundant coordinates automatically.
     
       The coordinates generated from that step will be generated from the current cartesian coordinates plus the constrained that you have set.
    
    
• to do a constrained geometry optimization, just run jobex [-ri -c 100 ...] in each of the subdirectories
• To collect the results, just grep for the value of the internal coordinates and the energy after the geometry optimizations have been performed

Hints:

• If you do not copy the initial input files to all directories, but write a script that starts with the first geometry optimization, and continues with the next step by using the coordinates of the step before, the scan might converge faster if the structure changes significantly during optimization.
• The 'trick' to define just one fixed internal coordinate and let ired find the correct new structure can also be used for other purposes. E.g., if you are looking for conformers and want to change just a dihedral angle, simply call define, fix the torsion, change its value to whatever you want, call ired and you will get the coordinates of the twisted structure.
• You are not restricted to one fixed internal coordinate. Define as many as you need and keep them fixed during the hole procedure. If you wrote a script that does a scan half-automatically, you can fix two internal coordinates, do a scan as described above for the first one, enter each of the directories, and run the same script for the second fixed internal coordinate. This will result in a 2D-scan.
• Instead of using the default module for geometry optimizations relax, you can use statpt by adding the option -statpt to the jobex call. statpt is often more stable and sometimes faster than the default.