ioutil.c

// File: ioutil.c
/*
 * Input/Output routines for pjayci.x
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "ioutil.h"
#include <mpi.h>

/* ================================================================= */
/* Fortran subroutines */
extern void _readmoints_(double *moints1, double *moints2,
			 long long int itype,
			 long long int orbitals,
			 long long int m1len,
			 long long int m2len,
			 double       *energy);

extern void _readnamelist_(long long int  nmlist,
			   unsigned char *nmlstr,
                           long long int *err);

extern void _readmocoef_(double *c, long long int clen);
/* ================================================================= */


/* check_for_file
 * --------------
 * Check if file exists and can be opened for desired operation.
 */
int check_for_file(char *filename, char *fileoperation)
{
        int iexists = 0;
        FILE *fptr = NULL;
        char error_message[80];
        fptr = fopen(filename, fileoperation);
        if (fptr == NULL){
                sprintf(error_message,"Cannot open file: '%s'; mode: '%s'!",
                       filename, fileoperation);
                iexists = 1;
                printf("Error: %s\n", error_message);
                return iexists;
        }
        fclose(fptr);
        return iexists;
}

/* checkinputfiles
 * ---------------
 * Check for the following input files:
 *  (1) = jayci.in    = input file containing &general and &diagalg namelists
 *  (2) = input.jayci = input file generated by jayci_exp.x
 *  (3) = moints      = molecular integral file
 *  (4) = det.list    = determinant input list
 * Returns integer value of missing input file.
 */
int checkinputfiles()
{
     /* local scalars
      * err = error handling */
     int err = 0;

     err = check_for_file("jayci.in","r");
     if (err != 0) return 1;
     err = check_for_file("input.jayci","r");
     if (err != 0) return 2;
     err = check_for_file("moints","r");
     if (err != 0) return 3;
     err = check_for_file("det.list","r");
     if (err != 0) return 4;

     return err;
     
}

/* find_str_count_in_file: find occurances of string in file.
 * -------------------------------------------------------------------
 * Input:
 *  string = string to search for
 *  stream = file stream
 */
int find_str_count_in_file(char *string, FILE *fptr)
{
        int result = 0; /* Count */
        char tmp[MAX_LINE_SIZE];

        rewind(fptr);
        while (fgets(tmp, MAX_LINE_SIZE, fptr) != NULL) {
                if ((strstr(tmp, string)) != NULL) {
                        result++;
                }
        }
        if (result == 0) {
                result = -1;
        }
        return result;
}

/* find_str_line: find and return position in FILE *stream where string
 * first occurs.
 * -------------------------------------------------------------------
 * Input:
 *  string = string to search for
 *  stream = file stream
 */
FILE *find_str_line(char *string, FILE *fptr)
{
        char tmp[MAX_LINE_SIZE];
        fpos_t ptr;

        rewind(fptr);
        fgetpos(fptr, &ptr);
        while (fgets(tmp, MAX_LINE_SIZE, fptr) != NULL) {
                if ((strstr(tmp, string)) != NULL) {
                        fsetpos(fptr, &ptr);
                        return fptr;
                }
                fgetpos(fptr, &ptr);
        }
        return fptr;
}

/*
 * print_array_2d: print 2d array
 */
void print_array_2d(double **array, int rows, int cols)
{
        int i, j;
        for (i = 0; i < rows; i++) {
                for (j = 0; j < cols; j++) {
                        fprintf(stdout, " %15.8lf ", array[j][i]);
                }
                fprintf(stdout,"\n");
        }
        return;
}

/*
 * print_wavefunction_info: print wavefunction information.
 */
void print_wavefunction_info(char *wfname, int nelecs, int norbs, int nfrzc,
                             int ndocc, int nactv, int nfrzv, int xlvl,
                             int nstates)
{
        printf("%s: electrons = %d orbitals = %d\n", wfname, nelecs, norbs);
        printf(" (Frozen core, DOCC, ACTV, Frozen virtual) = %d %d %d %d\n",
               nfrzc, ndocc, nactv, nfrzv);
        printf(" Number of states: %d\n", nstates);
        return;
}

/*
 * readdysoninput: read &dysonorbital namelist.
 * Output:
 *  states0 = states of anion to compute dyson orbitals
 *  states1 = states of neutral to compute dyson orbitals
 *  error   = error flag
 */
void readdysoninput(int *states0, int *states1, int maxst, int *nst0,
                    int *nst1, int *error)
{
        long long int dysonnml = 5; /* &dysonorbital namelist flag */
        char nmlstr[MAX_NAMELIST_SIZE][MAX_LINE_SIZE] = {{""},{""}};
        int i, j;
        *error = 0;
        /* Read namelist */
        readnamelist_(&dysonnml, nmlstr, &error);
        if (error != 0) return;
        /* Stream the input into the proper variables */
        for (i = 0; i < maxst; i++) {
                j = i + maxst;
                sscanf(nmlstr[i], "%d", &states0[i]);
                sscanf(nmlstr[j], "%d", &states1[i]);
                if (states0[i] != 0) *nst0 = (i + 1);
                if (states1[i] != 0) *nst1 = (i + 1);
        }
        return;
}
        
/* readwf0input: read wavefunction input for anion (0)
 * -------------------------------------------------------------------
 * Calls readnamelist which returns a character array
 *  nmlstr[0] = elec
 *  nmlstr[1] = orbs
 *  nmlist[2] = nfrozen
 *  nmlist[3] = ndocc
 *  nmlist[4] = nactive
 *  nmlist[5] = xlevel
 *  nmlist[6] = nfrzvirt
 *  nmlist[7] = nstates
 *
 * Output:
 *  elec = number of electrons in system (alpha + beta)
 *  orbs = number of orbitals in system (including forzen core)
 *  nfrozen = number of frozen core orbitals
 *  ndocc = number of doubly-occupied orbitals
 *  nactive = number of active orbitals
 *  xlevel = excitaion level (Default is 2)
 *  nfrzvirt = number of frozen virtual orbitals
 *  nstates = number of states
 *  err = error handling: n = missing variable n */
void readwf0input(int *elec,     int *orbs,   int *nfrozen,  int *ndocc,
	          int *nactive,  int *xlevel, int *nfrzvirt,
                  int *nstates,  int *err)
{
     /* local scalars
      * dywf0nml = namelist to read in */
     long long int dywf0nml=3;
     
     /* local arrays
      * nmlstr = namelist character arrays */
     char nmlstr[MAX_NAMELIST_SIZE][MAX_LINE_SIZE] = {{""},{""}};
     

     *err = 0;

     /* read namelist 1 */
     readnamelist_(&dywf0nml, nmlstr, &err);
     if (err != 0) return;

     /* stream the input into the proper variables */
     sscanf(nmlstr[0], "%d", elec);
     sscanf(nmlstr[1], "%d", orbs);
     sscanf(nmlstr[2], "%d", nfrozen);
     sscanf(nmlstr[3], "%d", ndocc);
     sscanf(nmlstr[4], "%d", nactive);
     sscanf(nmlstr[5], "%d", xlevel);
     sscanf(nmlstr[6], "%d", nfrzvirt);
     sscanf(nmlstr[7], "%d", nstates);

     return;
     
}

/* readwf1input: read wavefunction input for neutral (1)
 * -------------------------------------------------------------------
 * Calls readnamelist which returns a character array
 *  nmlstr[0] = elec
 *  nmlstr[1] = orbs
 *  nmlist[2] = nfrozen
 *  nmlist[3] = ndocc
 *  nmlist[4] = nactive
 *  nmlist[5] = xlevel
 *  nmlist[6] = nfrzvirt
 *  nmlist[7] = nstates
 *
 * Output:
 *  elec = number of electrons in system (alpha + beta)
 *  orbs = number of orbitals in system (including forzen core)
 *  nfrozen = number of frozen core orbitals
 *  ndocc = number of doubly-occupied orbitals
 *  nactive = number of active orbitals
 *  xlevel = excitaion level (Default is 2)
 *  nfrzvirt = number of frozen virtual orbitals
 *  nstates = number of states to read
 *  err = error handling: n = missing variable n */
void readwf1input(int *elec,     int *orbs,   int *nfrozen,  int *ndocc,
	          int *nactive,  int *xlevel, int *nfrzvirt,
                  int *nstates,  int *err)
{
     /* local scalars
      * dywf0nml = namelist to read in */
     long long int dywf1nml=4;
     
     /* local arrays
      * nmlstr = namelist character arrays */
     char nmlstr[MAX_NAMELIST_SIZE][MAX_LINE_SIZE] = {{""},{""}};
     

     *err = 0;

     /* read namelist 1 */
     readnamelist_(&dywf1nml, nmlstr, &err);
     if (err != 0) return;

     /* stream the input into the proper variables */
     sscanf(nmlstr[0], "%d", elec);
     sscanf(nmlstr[1], "%d", orbs);
     sscanf(nmlstr[2], "%d", nfrozen);
     sscanf(nmlstr[3], "%d", ndocc);
     sscanf(nmlstr[4], "%d", nactive);
     sscanf(nmlstr[5], "%d", xlevel);
     sscanf(nmlstr[6], "%d", nfrzvirt);
     sscanf(nmlstr[7], "%d", nstates);

     return;
     
}

/* readdaiinput: read diagonalization algorithm input.
 * -------------------------------------------------------------------
 * Calls readnamelist which returns a chracter array
 *  nmlstr[0] = maxiter
 *  nmlstr[1] = krymin
 *  nmlstr[2] = krymax
 *  nmlstr[3] = nroots
 *  nmlstr[4] = prediagr
 *  nmlstr[5] = refdim
 *  nmlstr[6] = restol
 *  nmlstr[7] = ga_buflen
 *
 * Output:
 *  maxiter = maximum iterations of davidson algorithm
 *  krymin  = minimum dimension of krylov space
 *  krymax  = maximum dimension of krylov space
 *  nroots  = number of roots to find
 *  prediagr= prediagonalization subroutine choice
 *  refdim  = intitial reference-space dimension (prediagonalization)
 *  restol  = convergence tolerance of residual
 *  buflen  = length of buffer in GA read of Hv=c
 *  err     = error handling: n = missing variable n */
void readdaiinput(int *maxiter,  int *krymin, int *krymax, int *nroots,
		  int *prediagr, int *refdim, double *restol, int *buflen,
                  int *err)
{
    /* .. local scalars ..
     * gnml = namelist to read in */
    long long int gnml = 2;

    /* .. local arrays ..
     * nmlstr = namelist character arrays */
    char nmlstr[MAX_NAMELIST_SIZE][MAX_LINE_SIZE] = {{""},{""}};

    *err = 0;

    readnamelist_(&gnml, nmlstr, &err);
    if (err != 0) return;

    /* stream input into proper values */
    sscanf(nmlstr[0], "%d",  maxiter);
    sscanf(nmlstr[1], "%d",   krymin);
    sscanf(nmlstr[2], "%d",   krymax);
    sscanf(nmlstr[3], "%d",   nroots);
    sscanf(nmlstr[4], "%d", prediagr);
    sscanf(nmlstr[5], "%lf",  restol);
    sscanf(nmlstr[6], "%d",   refdim);
    sscanf(nmlstr[7], "%d",   buflen);
    
    return;
}

/* readgeninput: read general wavefunction input.
 * -------------------------------------------------------------------
 * Calls readnamelist which returns a character array
 *  nmlstr[0] = elec
 *  nmlstr[1] = orbs
 *  nmlist[2] = nfrozen
 *  nmlist[3] = ndocc
 *  nmlist[4] = nactive
 *  nmlist[5] = xlevel
 *  nmlist[6] = nfrzvirt
 *  nmlist[7] = printlvl
 *  nmlist[8] = printwvf
 *
 * Output:
 *  elec = number of electrons in system (alpha + beta)
 *  orbs = number of orbitals in system (including forzen core)
 *  nfrozen = number of frozen core orbitals
 *  ndocc = number of doubly-occupied orbitals
 *  nactive = number of active orbitals
 *  xlevel = excitaion level (Default is 2)
 *  nfrzvirt = number of frozen virtual orbitals
 *  printlvl = print level
 *  printwvf = print wavefunctions (0: no; 1: yes)
 *  err = error handling: n = missing variable n */
void readgeninput(int *elec,     int *orbs,   int *nfrozen,  int *ndocc,
	          int *nactive,  int *xlevel, int *nfrzvirt, int *printlvl,
                  int *printwvf, int *err)
{
     /* local scalars
      * gnml = namelist to read in */
     long long int gnml=1;
     
     /* local arrays
      * nmlstr = namelist character arrays */
     char nmlstr[MAX_NAMELIST_SIZE][MAX_LINE_SIZE] = {{""},{""}};
     

     *err = 0;

     /* read namelist 1 */
     readnamelist_(&gnml, nmlstr, &err);
     if (err != 0) return;

     /* stream the input into the proper variables */
     sscanf(nmlstr[0], "%d", elec);
     sscanf(nmlstr[1], "%d", orbs);
     sscanf(nmlstr[2], "%d", nfrozen);
     sscanf(nmlstr[3], "%d", ndocc);
     sscanf(nmlstr[4], "%d", nactive);
     sscanf(nmlstr[5], "%d", xlevel);
     sscanf(nmlstr[6], "%d", nfrzvirt);
     sscanf(nmlstr[7], "%d", printlvl);
     sscanf(nmlstr[8], "%d", printwvf);

     return;
     
}

/* readmointegrals: Subroutine to read 1 and 2 electron integrals.
 * -------------------------------------------------------------------
 * Calls fortran subroutine readmoints()
 *
 * Input:
 *  itype     = type of integrals to read
 *  orbitals = MO's in system
 *  mofile   = name of molecular orbital file
 *  m1len    = length of moints1
 *  m2len    = length of moints2
 * Output:
 *  moints1  = 1-e integrals
 *  moints2  = 2-e integrals
 *  nuc_rep  = nuclear repulsion energy
 *  fcenergy = frozen-core energy */
void readmointegrals(double *moints1, double *moints2, int itype,
		     int orbitals, char *restrict moflname, int m1len,
		     int m2len, double *nuc_rep, double *fcenergy)
{
     long long int itype8, orbitals8, m1len8, m2len8;
     double energy[2];

     itype8 = (long long int) itype;
     orbitals8 = (long long int) orbitals;
     m1len8 = (long long int) m1len;
     m2len8 = (long long int) m2len;
     
     fprintf(stdout, "Calling readmoints_\n");
     fprintf(stdout, " Molecular integral file: %s\n", moflname);
     fprintf(stdout, " Type of integrals: %lld\n", itype8);
     fprintf(stdout, " 1-e integrals: %lld\n", m1len8);
     fprintf(stdout, " 2-e integrals: %lld\n", m2len8);
     /* call fortran subroutine */
     readmoints_(moints1, moints2, &itype8, &orbitals8, &m1len8,
		&m2len8, energy);
     *nuc_rep =  energy[0];
     *fcenergy = energy[1];
     return;
}

/*
 * substring: gets substring from string and returns pointer to said
 * substring.
 */
char *substring(char *string, int position, int length)
{
        char *pointer;
        int c;

        pointer = malloc(length+1);
        if (pointer == NULL) {
                printf("Unable to allocate memory.\n");
                exit(1);
        }
        for (c = 0 ; c < length ; c++) {
                *(pointer+c) = *(string+position);
                string++;
        }
        *(pointer+c) = '\0';
        return pointer;
}