Programming Your Own Input Routine

VIEWMOL can be easily extended to read outputs of other programmes. All you have to do is to write a subroutine which extracts the data from the output file. You have to use the routines for dynamical memory allocation. These are three routines:

• getmem(pointer,number_of_elements,size_of_element)
  This subroutine allocates number_of_elements times size_of_ele-
ment bytes memory from the operating system. pointer is a pointer on the allocated memory. The viewmol.par file contains a lot of pointer definitions for all global arrays. Use these pointers in calling getmem. After allocating memory with getmem you can use the arrays in the normal way (except if the array has more than one dimension).
• expmem(pointer,number_of_elements,size_of_element)
  This subroutine extends the memory previously allocated to number_
of_elements times size_of_elements bytes. Specifying less memory than in getmem reduces the size of the allocated memory. The contents of the memory which is not influenced by expanding/shrinking is not changed. You can use this subroutine when you find out while reading data that your arrays are too small. Then you call expmem with suitable parameters and continue with reading data.
• fremem(pointer)
  This subroutine frees memory previous allocated. Don't use this subroutine for the global arrays !

cphist copies the data from the three dimensional hist array (see below) to another destination or vice versa. This routine should be used when extending the memory for the optimization history data. You allocate first memory for a temporary array using getmem. Then you copy the contents of the array hist to the temporary array using cphist. Now you can expand the array hist with the routine expmem. Finally you copy back the contents of the hist array using cphist and free the memory of the temporary array with fremem.

cpmo copies the data from the MO array to another destination or vice versa. It works exactly the same way as cphist, but for two-dimensional arrays.

The following code fragment contains an example for using the dynamical memory allocation. You can also look in the data reading subroutines always provided.

c
c      Cartesian coordinates are saved in arrays x(), y() and
c      z(). The pointers to these arrays are ipx, ipy and ipz.
c      These arrays are real*4 arrays.
c
       mna=50                   ! default number of atoms
       call getmem(ipx,mna,4)   ! allocate memory for arrays
       call getmem(ipy,mna,4)
       call getmem(ipz,mna,4)
       i=1
    10 read(7,*) x(i),y(i),z(i) ! read data from file
       i=i+1
       if (i.gt.mna) then       ! more atoms than we have
                                ! expected
         mna=mna+50             ! make maximum number of
                                ! atoms larger
         call expmem(ipx,mna,4) ! and expand arrays
         call expmem(ipy,mna,4)
         call expmem(ipz,mna,4)
       endif
       goto 10                  ! continue with next atom

The second example demonstrates the use of the indarr and the
cphist routines.

c
c      Optimization history data are saved in three dimensio-
c      nal array hist(6,na,mnh). This array is a real*4 array.
c      Na is the known number of atoms, mnh is the still
c      unknown number of optimization history sets (ohs).
c
       pointer(iphelp,help)       ! define iphelp as pointer
                                  ! on array help
       dimension help(*)          ! make help being an array
       mnh=20                     ! default number of ohs
       call getmem(iphist,6*na*mnh,4) ! allocate memory for
                                      ! array hist
       i=1
    10 do 20 j=1,na
    20 read(7,*) (hist(indarr(k,j,i,6,na)),k=1,6)
                                  ! read data from file
       i=i+1
       if (i.gt.mnh) then         ! more ohs than we have
                                  ! expected
         call getmem(iphelp,6*na*mnh,4) ! allocate memory for
                                        ! temporary array
         call cphist(hist,help,mnh,mnh,na) ! copy hist to
                                           ! help
         mnh=mnh+20               ! make mnh larger
         call expmem(iphist,6*na*mnh,4) ! expand array hist
         call cphist(help,hist,mnh-20,mnh,na) ! copy help to
                                  ! hist with new dimensions
         call fremem(iphelp)      ! free array help
       endif
       goto 10                    ! next ohs

The data you have to read in in your own input routine are (all arrays are one dimensional, exceptions are noted):

• the cartesian coordinates
  They must be stored in the arrays x(), y() and z(). The array symbol() (character*2) must contain the atomic symbol. The variable radfac must be set to a value so that radfac*x(i) etc. gives the coordinate in Ångstrøm. The variable na  has to contain the number of atoms read in. All these data are essential.
• the title
  The character*80 variable title has to contain the title of the calculation if there is any.
• the wave numbers, intensities and normal coordinates
  The wave numbers in cm must be read into array wn(). The intensities must be read into array wi(). After reading all intensities one has to call subroutine relint which calculates relative intensities from the data in wi() so that the units of the intensities do not matter. If the symmetries of the normal vibrations can be found in the file they have to be read into array symmod() (character*4). If this information is not available the array should be set to a default (e. g. A1 or ?). The arrays ir() and raman() (logical) have to be set to .true. or .false. depending on whether the corresponding mode is IR or Raman active. If this information is not available both arrays should be set to .true. The variable ncart must contain the number of wave numbers read in.

  The normal coordinates must be read into array anc() (two dimensional) where the first index runs over the wave numbers and the second one over the 3*na cartesian coordinates of the atoms. After reading all data successfully the logical variable specex has to be set to .true. If these data could not found specex has to be .false.

• optimization history
  If the output contains informations about the history of a geometry optimization these data can be read into arrays energy(), gnorm() and hist (three dimensional). energy() and gnorm() must contain the energy and the gradient norm of each step of the geometry optimization. The array hist() must contain the coordinates and the gradients of each step of the optimization. The first index of this array runs from 1 to 6, where 1, 2 and 3 are the x, y and z component of the coordinates and 4, 5 and 6 the x, y and z component of the gradient. The second index counts the atoms and the third one the step of the geometry optimization. The number of steps from the geometry optimization has to be put into the variable nhist. If the geometry optimization converged the logical variable statio has to be .true. If there are gradients for the last step of the geometry optimization igrad has to be equal to nhist otherwise equal to -1.
• energy level diagram
  If the output contains informations about MO energies these data can be read into arrays eorb() and symorb(). eorb() must contain the energies in Hartree sorted so that eorb(1) contains the lowest and eorb(nbasfu) the highest MO energy. symorb() has to contain the symmetry labels for each MO. After reading the data the following logical variables have to be set to .true.: exmoe if orbital energies have been read in and exmos if orbital symmetries have been read in. The integer variable nbasfu has to contain the number of MO's.
• wave function related data
  For the drawing of basis functions, MO's, and electron densities you need basis functions, MO coefficients and occupation numbers. Since the handling of basis functions can be rather complicated, a function is supplied which stores the basis functions into the appropriate arrays. This function is called addatm and handles the basis functions per atom. Call it the following way with your basis set: call addatm(na,i,ishell,maxshe,maxset,inshell,alfas,coeff, exbas) na is the total number of atoms (as set in the coordinate input, see p. ), i is the number of the actual atom, ishell specifies how many contracted Gaussians atom i has, maxshe and maxset are the dimensions of the array inshell, inshell is an integer array of the dimensions (maxshe, maxset+2) which specifies the contraction scheme of the basis set (vide infra), alfas is an one-dimensional floating point array which contains the basis function exponents, coeff is a two-dimensional floating point array which contains the contraction coefficients (have to be 1. if only one primitive per contracted basis function), and exbas is a logical variable which will be set by addatm to .true. if storing the basis functions was successful. The contraction scheme is given in inshell in the following way: inshell(j,1) specifies how many sets of contraction coefficients per exponent are in the coeff array (usually 1, but with Pople basis sets there may be 2 - one for s and one for p). inshell(j,2) specifies how many primitives shell j has. inshell(j,3 ...) specify the angular quantum number for each set of coefficients (s=0, p=1, d=2, f=3, g=4). coeff has to be of the dimensions (maxset,*). If the basis set uses cartesian basis functions (6d, 10f, 15g) the logical variables pured, puref, and pureg have to be set to .false., with "true" basis functions (5d, 7f, 9g) these variables have to be set to .true.

  The MO coefficients have to be read into a two-dimensional floating point array cmo column by column. For the handling of dynamical memory allocation the function cpmo is provided which works exactly as cphist. If the MO coefficients have been successfully read in, the logical variable exmoc has to be set to .true. The MO coefficients have to be stored with respect to cartesian basis functions regardless whether these were used in the ab initio calculation. Most programmes will print the MO coefficients with respect to cartesian basis functions. If your programme does not, you have to figure out how your output can be transformed to cartesian basis functions. You may take a look into the function modcmo which does this transformation for TURBOMOLE and GAUSSIAN. This function also lists in which order your input function must store the MO coefficients for the different magnetic quantum numbers.

  The occupation number for each MO has to be read into the integer array nocc.

  The routines which calculate MO's and electron densities don't use symmetry. You can use them only if you have all MO coefficients for all atoms.

Your data reading subroutine has to be called from subroutine input (see source code). You have to specify a letter for your data format (currently d, g, m, o, p, and t are used). Then you can extend the input handler. The existence of your input file (the name you typed on the command line is in variable file) is checked before the corresponding input subroutine is called. You should open and close your file in the corresponding subroutine. If your file contains an unique identifier you may check it before calling your subroutine to make sure that you call the right one.