Geometry section of COLUMBUS The geometry optimization is performed by the GDIIS algorithm developed by Csaszar and Pulay (J. Mol. Structure, 114, 31 (1984) using a modified DIIS matrix. The geometry optimization is carried out in internal coordinates. The definition of these coordinates can be generated with the program INTC developed by Fogarasi et al. (G. Fogarasi, X.Zhou, P.W. Taylor and P. Pulay, JACS 114,8191 (1992). These definitions are stored in the file INTCFL and are not changed during the geometry optimization cycles. For a documentation of the INTCFL file see below. Each time GDIIS is executed one step of a geometry optimization cycle is performed. Input to GDIIS is the current geometry and cartesian gradient and output is a new geometry in the course of a energy minimum, saddle point or conical intersection search. This new geometry can be used for a new gradient calculation in the next optimization step. Documentation of the program G D I I S (Version: 20-Jul-03) 0. Introduction This program is designed to drive geometry optimization (search for equilibrium geometry, saddle point or conical intersection). The outer 'shell' of the program (communication with the quantum chemistry program) was designed and programmed by Peter G. Szalay and modified later by a few people. The inner 'shell', the optimization algorithm was taken from Peter Pulay's TX90 program system (version 1995. summer). BFGS option included by Peter Hoechtl, University of Vienna, April 1996. THIS PROGRAM CANNOT BE USED WITHOUT THE PERMISSION OF THE AUTHORS. 1. Motivation The GDIIS algorithm has proven to be one of the most efficient geometry optimization algorithms. The aim of this version was to utilize the advanced algorithm present in TX90 for other quantum chemistry program packages. Therefore the input needed for GDIIS is designed to be as general as possible to allow easy porting to any program systems. At the same time, we tried to avoid the necessity of changing GDIIS in the porting process. Two main files are necessary for communication: gdiisin includes input lines communicating geometry, energy, forces, internal coordinates and control parameters (including file names) to GDIIS gdiismg this file is used to pass information like new geometry, convergence information from GDIIS back to the quantum chemistry program (the name of this file might be changed as an option in gdiisin file, see below) By the curse of time this version became an 'integrated' part of COLUMBUS and it makes use of several files for communication with other programs of COLUMBUS. 2. The gdiisin file The records of the gdiisin file start with a keyword telling the program what information is included on the current record (eventually the proceeding records). The order of the keywords is arbitrary except when noted. It is possible to have information such as geometry, energy, forces on separate file(s). In this case the corresponding control record contains the file name. If the file name is gdiisin (default) the information will be read from the next record of the gdiisin file. 2.1. Keyword ener Record format: ener enerfl where ener : keyword (character*4) enerfl : file containing the energy of the current iteration (character*10, default: convfl) Purpose: to pass total energy to GDIIS Notes: -the keyword is required if 'sadd' option is used 2.2 Keyword geom Record format: geom na geomfl geomfm geomst (a4,i4,a10,a30,a30) where geom : keyword (character*4) na : number of atoms (integer*4). 0 is sufficient if the geometry is written till the end of file geomfl. geomfl : file containing the current geometry (character*10, default: gdiisin) geomfm : format of the records of geomfl file (character*30, default: '(2x,A2,6x,4F10.6)') geomst : style of the geomfl file. The available styles are listed below (character*30, default: texas) Next records in geomfl file if geomst='texas' or 'TEXAS' for i=1,na (if na=0, until the end of file is reached): symb(i),xiat,xa(1,i),xa(2,i),xa(3,i) where symb(i): symbol for atom i xiat: periodic number of atom i xa(1-3,i): x,y,z coordinate of atom i else: no other style is available presently (new options can be added by editing the subroutine readgeom) Purpose: to pass cartesian coordinates to GDIIS Notes: -the keyword is required 2.3 Keyword forc (only after keyword geom!) Record format: forc forcfl forcfm forcst (a4,1x,a10,a30,a30) where forc : keyword (character*4) forcfl : file containing the current cartesian forces (character*10, default: gdiisin) forcfm : format of the records of forcfl file (character*30, default: '(3F10.6)' forcst : style of the forcfl file. The available styles are listed below (character*30, default: texas) Next records in forcfl file if forcst='texas' or 'TEXAS' f(i) for i=1,3*na where f(3j-2),f(3j-1),f(3j): the x,y,z component of the forces on atom j IMPORTANT: forces (i.e. negative gradients) are required!!!! else: no other style is available now Purpose: to pass cartesian forces to GDIIS Notes: -the keyword is required 2.4 Keyword intc Record format: intc nc intcfl (a4,i4,a10) where intc : keyword (character*4) nc : number of internal coordinates. Nonzero value means that internal coordinates are written on the next nc records (default: 0) intcfl : nc=0 file intcfl contains the definition of the internal coordinates nc.ne.0 internal coordinates are on the lines following directly the key word line Next nc records are read with format A80 Purpose: to pass the internal coordinates and its file names. The internal coordinates are in the style of TX90 and described also below Notes: -the keyword is not required if the internal coordinates are on file intcfl or generated by keyword 'auto' 2.5 Keyword fmfl Record format: fmfl nc fmatfl (a4,i4,a10) where fmfl : keyword (character*4) nc : number of internal coordinates. Nonzero value means that force constants are written on the next nc records (default: 0) fmatfl : file containing (nc=0) or reserved (n.ne.0) for force constants. Next nc records are read with format A80 (if nc.ne.0) Purpose: to pass the force constants and its file names Notes: -the keyword is optional 2.6 Keyword form NOT SUPPORTED IN THE PRESENT VERSION Purpose: to pass information about previous iterations Notes: -the keyword is optional 2.7 Keyword prin Record format: prin where prin : keyword (character*4) Purpose: activates higher level of printing Notes: -the keyword is optional 2.8 Keyword fixc Record format: fixc x (a4,i4) where fixc : keyword (character*4) x : number of the internal coordinate to be fixed Purpose: fixes internal coordinate x in the optimization procedure Notes: - the keyword can be used several times - the keyword is optional 2.9 Keyword fint NOT SUPPORTED IN THE PRESENT VERSION Record format: fint where fint : keyword (character*4) Notes: -the keyword is not yet installed 2.10 Keyword auto NOT SUPPORTED IN THE PRESENT VERSION Record format: auto where auto : keyword (character*4) Purpose: invoke automatic procedure to generate internal coordinates by intc.x Notes: - the keyword can be used only the first time GDIIS is called!!! - the keyword is not supported in the present version 2.11 Keyword sadd Record format: sadd where sadd : keyword (character*4) Purpose: searching of first order saddle points (for details see documentation in TX90) Notes: - the keyword is optional, the bfgs option is switched off automatically 2.12 Keyword shft Record format: shft shftpt where shft : keyword (character*4) shftpt: shift point for torsions (default: 1.57079632795=pi/2) Purpose: changes the shift points for torsions (for details see documentation in TX90) Notes: - the keyword is optional 2.13 Keyword maxi Record format: maxi coormx (a4,6x,f10.6) where maxi : keyword (character*4) coormx : maximum change permitted in the internal coordinates (default: 0.3) Purpose: sets the maximum change permitted in the internal coordinates Notes: - the keyword is optional 2.14 Keyword opti NOT USED IN THE PRESENT VERSION Record format: opti mxopti (a4,i4) where opti : keyword (character*4) mxopti : maximum number of iterations (default: 20) Purpose: sets the maximum number of iterations Notes: - the keyword is not used in the present version 2.15 Keyword maxc NOT USED IN PRESENT VERSION; SEE KEYWORDS: cmax, crms,fmax, frms Record format: maxc check1 chmax (a4,1x,a5,1x,f20.10) where maxc : keyword (character*4) check1 : type of convergence criteria (character*5, default: energ) possible values are: energ,coord,force chmax : convergence criteria(f20.10) Purpose: to set the type and size of the convergence criteria. Default type is an approximate energy lowering based convergence check. The other two options available are based on coordinate change or maximum internal force criterion. Notes: - the keyword is not used in the present version 2.16 Keyword gdii Record format: gdii mgdi (a4,i4) where gdii : keyword (character*4) mgdi : number of point required to start gdiis interpolation (default: 2) Purpose: indicates that the gdiis algorithm is used instead of the default simple relaxation and sets how many points are necessary to perform gdiis extrapolation the first time. Notes: - the keyword is optional 2.17 Keyword fdia Record format: fdia (a4) where fdia : keyword (character*4) Purpose: indicates that diagonal force constant matrix is used Notes: - fdia has no effect in this COLUMBUS version 2.18 Keyword bfgs Record format: bfgs (a4) where bfgs : keyword (character*4) Purpose: BFGS update of inverse Hessian. The inverse Hessian is stored in the file 'hessianinv'. The matrix is stored in square form. To be used for energy minimum search only. Notes: If the file 'hessianinv' exists it will be used. If 'hessianinv' does not exist the 'hessian' file is used, if it does not exist the diagonal elements of the hessian will be read from the input (this is the file 'intcfl' or a separate file defined by the key word fmfl) 2.19 Keyword punc Record format: punc (a4) where punc : keyword (character*4) Purpose: ? Notes: - not installed 2.20 Keyword angs Record format: angs (a4) where angs : keyword (character*4) Purpose: indicates that the coordinates are in Angstroms instead of Bohr. Notes: - optional - must proceed option 'geom' 2.21 Keyword nuse Record format: nuse itopti where iter : keyword (character*4) itopti : maximum number of geometries and gradients used in the gdiis interpolation Default value is 1, which means no diis at all (i4). Purpose: to control the inetrpolation procedure Notes: - 2.22 Keyword mess Record format: mess gdiismg (a4,1x,10a) where mess : keyword (character*4) gdiismg : file to be used for passing information back from GDIIS (default: gdiismg) Purpose: to set the file name of message file Notes: - optional 2.23 Keyword cond Record format: cond condamp (a4,6x,f10.6) where cond : keyword (character*4) condamp : (6x,f10.6) (default 0.0) Purpose: damping value to be used in subroutine congrad, a value .ne. zero (e.g. 1e-7) is necessary when the number of internal coordinates is less than 3N-6 (empirical observation by H. Lischka, Vienna) 2.24 Keyword lonl Record format: lonl (a4,i4) where lonl : keyword (character*4) lonl : (a4,i4) (default 0) Purpose: default=0; if = 1, preform only linear internal -> Cartesian transformation step (equal to 1 iteration) 2.25 Keyword lalw Record format: lalw (a4,i4) where lalw : keyword (character*4) lalw : (a4,i4) (default 0) Purpose: default=0; if = 1, if the conventional iterative procedure has not converged in maxiter steps, perform linear step 2.26 Keyword cmax Record format: (a4,6x,f20.10) cmax value where cmax : keyword (character*4) value : f20.10 (default 1.0d-03) Purpose: convergence control, maximal coordinate change threshold 2.27 Keyword crms Record format: (a4,6x,f20.10) crms value where crms : keyword (character*4) value : f20.10 (default 1.0d-03) Purpose: convergence control, maximal coordinate square norm threshold 2.28 Keyword fmax Record format: (a4,6x,f20.10) fmax value where fmax : keyword (character*4) value : f20.10 (default 1.0d-03) Purpose: convergence control, maximal force threshold 2.29 Keyword frms Record format: (a4,6x,f20.10) frms value where frms : keyword (character*4) value : f20.10 (default 1.0d-03) Purpose: convergence control, maximal force square norm threshold 2.30 Keyword for1 Record format and description: same as for keyword force (2.15) Purpose: to pass Cartesian forces of state 1 to GDIIS Notes: -the keyword is required for 'coni' and 'avoi' calculations 2.31 Keyword for2 Record format and description: same as for keyword force (2.15) Purpose: to pass Cartesian forces of state 2 to GDIIS Notes: -the keyword is required for 'coni' and 'avoi' calculations 2.32 Keyword coup Record format and description: same as for keyword force (2.15) Purpose: to pass Cartesian coupling vectors to GDIIS Notes: -the keyword is optional 2.33 Keyword ene1 Record format: ene1 enerfl (a4,1x,10a) where ener : keyword (character*4) enerfl : file containing the energy of the current iteration (character*30, default: convfl) Next record in enerfl file etot where etot: is the energy of the current iteration (real*8) Purpose: to pass total energy of state1 to GDIIS Notes: -does not work: energy will be read from file 'energy'!! 2.34 Keyword ene2 Record format: ene2 enerfl (a4,1x,10a) where ener : keyword (character*4) enerfl : file containing the energy of the current iteration (character*30, default: convfl2) Next record in enerfl file etot where etot: is the energy of the current iteration (real*8) Purpose: to pass total energy of state2 to GDIIS Notes: -does not work: energy will be read from file 'energy'!! 2.35 Keyword zeta Record format: zeta zetafile (a4,1x,10a) where zeta : keyword (character*4) zetafile: file containing the lagrangians (zeta1 and zeta2) iteration (character*30, default: gdiisin) Description of zetafile: zeta1,zeta2 (free format) if the file does not exists, zeta1=-1. and zeta2=0 are the default values Purpose: to pass zeta1 and zeta2 to the program Notes: - 2.36 Keyword coni Record format: coni (a4) where coni : keyword (character*4) Purpose: indicates that conical intersection search should be performed Notes: - the keyword is optional 2.36 Keyword avoi Record format: avoi (a4) where avoi : keyword (character*4) Purpose: indicates that avoided crossing search should be performed Notes: - the keyword is optional (NOT TESTED!!!) 3. The gdiismg file 3.1 first record: information line Record format: excitecode where excitecode: excite information from GDIIS (A80) possible values are: - 'grad required; new geom follows' - 'stop; maximum number of iteration reached' - 'converged; convergence reached' 3.2 Next na records (na is the number of atoms): if geomst='texas' or 'TEXAS' for i=1,na symb(i),xiat,xa(1,i),xa(2,i),xa(3,i) where symb(i): symbol for atom i xiat: periodic number of atom i xa(1-3,i): x,y,z coordinate of atom i (unit is the same as in gdiisin; see keyword angs) else: no other style is available presently 4. The intcfl file contains the definition of the internal coordinates and the diagonal force constants. 4.1. internal coordinates For each elementary internal coordinate a card is read in the following format: kw, scale, coeff, type, number, a, b, c, d (A4, F6.0,....) kw may be either the character K (in column 1) or blank. If it is K, this shows that the present coordinate begins a new, independent internal coordinate. If it is blank then the coordinate is interpreted as the continuation of a composite linear coordinate begun earlier. Any other character terminates the input. scale is a scale factor for the total coordinate (only if there is K in column 1). Blank and zero is interpreted as 1.0 . coef is the linear combination factor of each elementary internal coordinate. type coordinate type: stre, invr, bend, out, tors, lin1, lin2 a,b,c,d participating atoms stre stretch coordinate, order of a and b arbitrary invr inverse stretch?, order of a and b arbitrary bend bond angle deformation, a and b are end atoms, c is the apex atom out out-of-plane deformation, atom a out of the bcd plane, c is the central atom, coordinate is positive if a is on the same side of the plane as the vector product db*bc. b and c can be exchanged but this changes the sign of the coordinate. tors the coordinate is defined as the angle of the planes abc and bcd. Note that abcd and dcba are equivalent; the sign of the torsion does not depend on the direction of the numbering. Note that the value of the coordinate varies between -pi/2 to 3pi/2, not between -pi/2 to pi/2 lin1 colinear bending of the linear chain of atoms a-c-b in the plane which contains d. The sign is positive if a and b move towards d. lin2 is the perpendicular linear bending. It is positive if a ad b move in the direction of the vector product cd*ca. 4.2 diagonal force constants: Record format: (8f10.7) (Format) 5. Files used filename filenumber purpose gdiisin inp = 30 input file gdiisls iout = 8 output file gdiissm icond= 9 summary file geomfl ngeom= 71 input file contains cartesian coordinates forcfl nforc= 72 input file contains cartesian forces intcfl nintc= 73 input file for definition of internal coordinates fmatfl nfmat= 74 input file contains force constants formfl nform= 75 enerfl nener= 76 input file contains calculated energy, and convergence informations gdiisfl inpf3= 77 output file stores internal coord. of prev. iterations gdiismg nmess= 78 output file for communication with other programs hessian nhmat= 79 output file stores formatted hessian hessianinv nhessinv=80 input file contains formatted inv hessian curr_iter ncuriter=81 contains the current iteration number 6. Current implementations: This version works only with the COLUMBUS program package. (There is a separate version for ACES2) 7. Final remarks: Any comment on this document or the program should be addressed to Peter G. Szalay (szalay@para.elte.hu) Concerning COLUMBUS implementations contact Hans Lischka (lischka@itc.univie.ac.at)