diff --git a/src/basis/GNUmakefile b/src/basis/GNUmakefile new file mode 100644 index 0000000..b5fb2df --- /dev/null +++ b/src/basis/GNUmakefile @@ -0,0 +1,15 @@ + HEADERS = bas.h basP.h bas_staticP.h + OBJ = basis.o newbasis.o bas_input.o + LIBRARY = libbasis.a + LIB_TARGETS = testbasis testbasis.o + +include ../config/makefile.h +include ../config/makelib.h + +testbasis: testbasis.o $(LIBRARY) + $(FC) $(FFLAGS) -o $@ testbasis.o $(LIBS) + + +basis.o: basP.h geobasmapP.h basdeclsP.h +basP.h: bas_staticP.h + @touch basP.h \ No newline at end of file diff --git a/src/basis/bas.h b/src/basis/bas.h new file mode 100644 index 0000000..7bb6cfd --- /dev/null +++ b/src/basis/bas.h @@ -0,0 +1,38 @@ +#ifndef _BAS_H +#define _BAS_H + +#include + +#ifdef __cplusplus +extern "C" { +#endif + +void bas_input(FILE *rtdb); +void bas_input_body(int basis, bool osegment); +bool gbs_map_clear(int basisin); + + bool bas_create(); + bool bas_destroy(); + bool bas_check_handle(int basisin, char *msg); + bool bas_321g_load(); + bool bas_print(); + bool bas_rtdb_load(); + bool bas_rtdb_store(); + bool bas_high_angular(); + bool gbs_map_print(); + bool bas_continfo(); + bool bas_numcont(); + bool bas_numbf(); + bool bas_get_exponent(); + bool bas_get_coeff(); + bool bas_set_exponent(); + bool bas_set_coeff(); + bool bas_print_all(); + bool bas_version(); + void bas_add_ucnt(int, char *, int, int, int, double*, double*, int); + +#ifdef __cplusplus +} +#endif + +#endif // _BAS_H_ diff --git a/src/basis/basP.h b/src/basis/basP.h new file mode 100644 index 0000000..196d8dd --- /dev/null +++ b/src/basis/basP.h @@ -0,0 +1,194 @@ +#ifndef _BASP_H +#define _BASP_H + +/* + basis set object/api + Rick A. Kendall and Robert J. Harrison (March 1994) + + What is the minimum basis set informaton? + + What is a basis set? A basis set is a set of tags nominally + associated with an atomic center through a geometry specification + or geometry object. What needs to be stored is the basis set for + the unique tags ("atoms"). The concept of shell and general + contraction confuses the issue somewhat but in the limit of a + segmented basis set the "shell" concept is the same as the + "general contraction." The basis set object and its interaction + with the integral API is predicated upon this assumption. + + A pseudo input deck with the minimum information is as follows. + Basis Set: Name (as on the rtdb) + ntags (number of unique tags for which basis set + information is supplied) + nucont (total number of unique contractions in basis) + nprim_t (total number of primitives in basis) + ncoef_t (total number of coeffs in basis) + foreach tag (ntags of them) + tag (character string identifier of tag) + number_of_contractions on tag + nprim_c in contractions on tag + ncoeff in contractions on tag + first contraction of tag + last contraction of tag + foreach contraction (number_of_contractions of them) + itype, nprim, ngen, iexptr, icoeff, tag_cont_is_on + foreach nprim in a contraction + ex(1), coeff(1,...) (ngen contractions) + + The tag ("atomic") information pseudo-data structure is as follows: + Integer num_cont ! Number of contractions on tag + Integer nprim_tag ! Number of primitive exponents on tag + Integer ncoeff_tag ! Number of primitive coeffs on tag + Integer ifirst_cont ! first contraction on tag + Integer ilast_cont ! last contraction on tag + Integer itype(num_cont) ! type of contraction 0=s,1=p,2=d... + ! later -1=sp,-2=spd etc. + Integer nprim(num_cont) ! number of primitives in each cont. on tag + Integer ngen(num_cont) ! number of general conts in each cont. + ! 1= segmented basis, >1 general cont. + Integer iexpt(num_cont) ! pointer into linearized real*8 array for + ! first exponent + Integer icoeffpt(num_cont) ! pointer into linearized real*8 array for + ! first coefficient + Integer itag(num_cont) ! tag identifier for contraction + ! (redundant for just atomic info) + double exndcf[nprim_tag+ncoeff_tag] ! linearized real*8 array of + ! exponents and coefficients + + The basis set information is the number of tags and the above tag + information. + + unique information only!!!!!!! + + The basis set information pseudo-data structure is as follows: + char bs_name[256] ! (as "mo basis") + int num_tags ! number of tags in basis + int num_cont_total ! number of conts in basis + int num_prim_total ! number of prims in basis + int num_coeff_total ! number of coefs in basis + char tags[num_tags][16] ! character string of tags + int num_cont[num_tags] ! Num of conts on tag + int nprim_tag[num_tags] ! Num of prim exponents on tag + int ncoeff_tag[num_tags] ! Num of prim coeffs on tag + int itype[num_cont] ! type of contraction 0=s,1=p,2=d... + ! later -1=sp,-2=spd etc. + int nprim[num_cont] ! num of prims in each cont. on tag + int ngen[num_cont] ! num of general conts in each cont. + ! on each tag + ! 1= segmented basis, >1 general cont. + int iexpt[num_cont] ! pointer into linearized real*8 array + ! for first exponent of each cont. + ! on each tag + int icoeffpt[num_cont] ! pointer into linearized real*8 + ! array for first coefficient of + ! each cont. on each tag + int itag[num_cont] ! tag identifier for contraction + ! (redundant for just atomic info) + double exndcf[] ! linearized real*8 array of + ! exponents and coefficients + Note: exndcf does not carry the num_cont label because it is a + linearized real*8 array and the added dimensionality is + handled by proper evaluation of the pointer arrays + iexpt and icoeffpt. + + The above data structure is too cumbersome to efficiently store to + and read from the run-time-data-base (one call per array). The + integer and real*8 data needs to be linearized with appropriate + informaton accessable by pointers arrays. + + The more appropriate "basis" data structure is as follows: + char bs_name[256] ! as "mo basis" + char tags[num_tags][16] ! character string of tags + int infbs_head[4] ! header information + ! 1 = num_tags ! num of tags in basis + ! 2 = num_cont_total ! num of conts in basis + ! 3 = num_prim_total ! num of prims(ex) in basis + ! 4 = num_coeff_total ! num of coeffs in basis + int infbs_tags[5][num_tags] + ! 1 = num_cont ! Num of conts on tag + ! 2 = nprim_tag ! Num of prim exponents on tag + ! 3 = ncoeff_tag ! Num of prim coeffs on tag + ! 4 = first cont ! first contaction on tag + ! 5 = last cont ! last contaction on tag + int infbs_cont[6][num_cont] + ! 1 = itype ! type of contraction 0=s,1=p,2=d... + ! later -1=sp,-2=spd etc. + ! 2 = nprim ! num of prims in each cont. on tag + ! 3 = ngen ! num of general conts in each cont. + ! on each tag + ! 1= segmented basis, >1 general cont. + ! 4 = iexpt ! pointer into linearized real*8 array + ! for first exponent of each cont. + ! on each tag + ! 5 = icoeffpt ! pointer into linearized real*8 + ! array for first coefficient of + ! each cont. on each tag + ! 6 = itag ! cunique tag number + double exndcf[] ! linearized real*8 array of + ! exponents and coefficients + + The above data structure now must handle multiple basis sets. + The "multiple basis set" data structure is as follows: + char bs_name[nbasis][256] ! as "mo basis" + char tags[num_tags][nbasis][16] ! character string of tags + int infbs_head[4][nbasis] + ! 1 = num_tags ! num of tags in basis + ! 2 = num_cont_total ! num of conts in basis + ! 3 = num_prim_total ! num of prims(ex) in basis + ! 4 = num_coeff_total ! num of coeffs in basis + int infbs_tags[5][num_tags][nbasis] + ! 1 = num_cont_tag ! Num of conts on tag + ! 2 = nprim_tag ! Num of prim exponents on tag + ! 3 = ncoeff_tag ! Num of prim coeffs on tag + ! 4 = first cont ! first contraction on tag + ! 5 = last cont ! last contraction on tag + int infbs_cont[6][num_cont_total][nbasis] + ! 1 = itype ! type of contraction 0=s,1=p,2=d... + ! later -1=sp,-2=spd etc. + ! 2 = nprim ! num of prims in each cont. on tag + ! 3 = ngen ! num of general conts in each cont. + ! on each tag + ! 1= segmented basis, >1 general cont. + ! 4 = iexpt ! pointer into linearized real*8 array + ! for first exponent of each cont. + ! on each tag + ! 5 = icoeffpt ! pointer into linearized real*8 + ! array for first coefficient of + ! each cont. on each tag + ! 6 = tag number ! unique center lexical index + double exndcf[][nbasis] ! linearized real*8 array of + ! exponents and coefficients + + actual names are set to protect name space between the geom and + basis objects +*/ +#define BASIS_HANDLE_OFFSET ((0-565)) + +// static dimension information for common blocks +#include "bas_staticP.h" + +// leading dimensions of compressed arrays +#define ndbs_tags 5 +#define ndbs_ucont 6 +#define ndbs_head 4 + +// stored structures +char bs_name[nbasis_bsmx][256]; +char bs_tags[ntags_bsmx][nbasis_bsmx][16]; +double exndcf[mxbs_exndcf][nbasis_bsmx]; +int infbs_head[ndbs_head][nbasis_bsmx]; +int infbs_tags[ndbs_tags][ntags_bsmx][nbasis_bsmx]; +int infbs_cont[ndbs_ucont][nucont_bsmx][nbasis_bsmx]; + +// in-core structures +double bsversion; +char bs_trans[nbasis_bsmx]; +char bs_names_rtdb[nbasis_rtdb_mx][256]; +int len_bs_name[nbasis_bsmx]; +int len_bs_trans[nbasis_bsmx]; +int len_bs_rtdb[nbasis_rtdb_mx]; +int nbasis_rtdb; +int angular_bs[nbasis_bsmx]; +bool bsactive[nbasis_bsmx]; + +#endif // _BASP_H diff --git a/src/basis/bas_input.c b/src/basis/bas_input.c new file mode 100644 index 0000000..4892d93 --- /dev/null +++ b/src/basis/bas_input.c @@ -0,0 +1,273 @@ +#include +#include + +#include "rtdb.h" +#include "context.h" +#include "geom.h" +#include "bas.h" +#include "inp.h" + +void bas_input(FILE *rtdb) { + /* + basis [] [library ] [file ] \ + [spherical|cartesian] [segment] [print] + + tag library [file ] + tag + + ... + end basis + + parse the main directive + */ + + int nopt = 6; + char opts[6][10] = {"spherical", "cartesian", "segment", "library", "file", "print"}; + char test[255], name[255], filename[255], standard[255]; + bool status, ospherical, osegment, oprint; + int ind, basis; + for (int i = 0; i < nopt; i++) { + opts[i][9] = '\0'; + } + + // Check is a basis directive and read in name of the basis + inp_set_field(0); + status = inp_a(test); + if ((!status) || (!inp_compare(false, test, "basis"))) { + goto L10000; + } + + // Parse rest of basis directive line + name[0] = '\0'; + filename[0] = '\0'; + test[0] = '\0'; + standard[0] = '\0'; + ospherical = false; + osegment = false; + oprint = false; + +L10: + if (inp_a(test)) { + if (!inp_match(nopt, false, test, opts, ind)) { + // Not a recognized option ... the name of the basis or an error + if ((name[0] != '\0') || (inp_cur_field() != 2)) { + printf(" bas_input: basis name must be first option\n"); + goto L10000; + } + strcpy(name, test); + goto L10; + } + switch (ind) { + case 0: // spherical + ospherical = true; + goto L10; + case 1: // cartesian + ospherical = false; + goto L10; + case 2: // segment + osegment = true; + goto L10; + case 3: // library + if (!inp_a(standard)) { + goto L10000; + } + goto L10; + case 4: // file + if (!inp_a(filename)) { + goto L10000; + } + goto L10; + case 5: // print + oprint = true; + goto L10; + } + } + + // Now check reality against input + if (standard[0] != '\0') { + printf(" Standard basis %s\n", standard); + errquit("bas_input: standard basis set not yet", 0); + } + if (ospherical) { + errquit("bas_input: spherical harmonics not yet", 0); + } + + // Open a new basis set to receive the new data + if (name[0] == '\0') { + strcpy(name, "mo basis"); + } + if (!bas_create(basis, name)) { + errquit("bas_input: failed to create basis", 0); + } + + // Here will soon process reading standard basis sets + + // Now left with reading in from the input additional specifications + // for basis functions or standard sets on specific tags + bas_input_body(basis, osegment); + + // Now have processed the entire basis directive. Print out + // info if desired, write it to the data base, tidy up and go home + if (oprint) { + if (!bas_print(basis)) { + errquit("bas_input: print failed", 0); + } + } + + if (!bas_rtdb_store(rtdb, name, basis)) { + errquit("bas_input: failed to store basis", 0); + } + + if (!bas_destroy(basis)) { + errquit("bas_input: bas_destroy failed", 0); + } + + return; + +L10000: + printf(" basis [] [library ] \\\n"); + printf(" [file ] [spherical|cartesian] [segment]\n"); + errquit("bas_input: invalid format for basis directive", 0); +} + +void bas_input_body(int basis, bool osegment) { +/* + +c Read the body of a basis directive that describes the +c tags/exponents/contraction coefficients +c +c +c tag library [file ] +c tag +c +c ... +c end basis +*/ + char tag[16], cont_type[16]; + int nltypes = 7, nsptypes = 2, nopts = 2; + int cont_max = 20, prim_max = 20; + double expnt[20], coeff[20][20]; + char ltypes[7] = {'s', 'p', 'd', 'f', 'g', 'h', 'i'}; + char sptypes[2] = {'sp', 'l'}; + int spvalues[2] = {-1, -1}; + char opts[2][8] = {"library", "file"}; + int spvalues[2] = {-1, -1}; + int l_value, ngen, iprim, nprim, i, ind; + +// Input a new line +L10: + if (!inp_read()) { + errquit("bas_input_body: premature EOF", 0); + } + +// Start parsing current line +L20: + inp_set_field(0); + if (!inp_a(tag)) { + goto L10000; + } + + if (inp_compare(false, "end", tag)) { + goto L9000; // End of basis directive + } + + if (!inp_a(cont_type)) { + goto L10000; + } + + if (inp_match(nltypes, false, cont_type, ltypes, ind)) { + // The contraction is a simple shell + l_value = ind - 1; + } else if (inp_match(nsptypes, false, cont_type, sptypes, ind)) { + // The contraction is an sp-type shell + l_value = spvalues[ind]; + } else if (inp_match(nopts, false, cont_type, opts, ind)) { + /* + It is actually an option to input a standard basis + + Don't bother parsing this yet + */ + errquit("bas_input_body: no standard basis sets yet", 0); + goto L10; // Process the next input line + } else { + // Only god and the user knows what was intended + goto L10000; + } + + // Fall thru to here to read in a set of contraction coefficients + if (!inp_read()) { + goto L10000; + } + + ngen = inp_n_field() - 1; + if (ngen < 1) { + goto L10000; + } + if (ngen > cont_max) { + errquit("bas_input_body: too many contractions - increase cont_max", cont_max); + } + + for (iprim = 0; iprim < prim_max; iprim++) { + if (!inp_f(expnt[iprim])) { + // If cannot read the first field as an exponent then + // it is the end of this contraction + + goto L30; + } else { + if ((inp_n_field() - 1) != ngen) { + printf(" bas_input_body: no. of coefficients?\n"); + goto L10000; + } + for (i = 0; i < ngen; i++) { + if (!inp_f(coeff[iprim][i])) { + printf(" bas_input_body: failed reading coefficient\n"); + goto L10000; + } + } + if (!inp_read()) { + goto L10000; + } + } + } + errquit("bas_input_body: too many primitives in contraction", prim_max); + +L30: + nprim = iprim - 1; + if (nprim <= 0) { + errquit("bas_input_body: no primitives?", nprim); + } + + // Now have tag, contraction type, no. of contractions, no. of prims, + // exponents, coeffs. Shove this lot into the basis set. + // bas_add_ucnt -> adds a new general contraction on the specified tag. + // If the tag is not present it will also add that. + if (osegment) { + // Add contractions one-at-a-time to force segmentation + for (i = 0; i < ngen; i++) { + if (!bas_add_ucnt(basis, tag, l_value, 1, nprim, expnt, coeff[0][i], prim_max)) { + errquit("bas_input_body: bas_add_ucnt failed!!", 0); + } + } + } else { + // Add as a single general contraction + if (!bas_add_ucnt(basis, tag, l_value, ngen, nprim, expnt, coeff, prim_max)) { + errquit("bas_input_body: bas_add_ucnt failed!!", 0); + } + } + + // Have already read in the next line ... parse it + goto L20; + +// Have read in all of the basis set info. +L9000: + return; + +L10000: + printf(" basis directive body format is:\n"); + printf(" tag library [file ]\n"); + printf(" tag \n"); + printf(" \n"); + printf(" ... \n"); + printf(" end basis\n"); + errquit("bas_input_body: format error in the input", 0); +} diff --git a/src/basis/bas_staticP.h b/src/basis/bas_staticP.h new file mode 100644 index 0000000..39052a6 --- /dev/null +++ b/src/basis/bas_staticP.h @@ -0,0 +1,14 @@ +#ifndef _BAS_STATICP_H +#define _BAS_STATICP_H + +// Maximum parameter definitions for static "in-core" data structure +const int nbasis_bsmx = 5; +const int nbasis_rtdb_mx = 10 * nbasis_bsmx; +const int ntags_bsmx = 10; +const int nucont_bsmx = 150; +const int mxbs_exndcf = ((300 + 500) * ntags_bsmx); + +const int nat_mx = 1000; +const int ncont_mx = (nucont_bsmx * 4); + +#endif // _BAS_STATICP_H diff --git a/src/basis/basdeclsP.h b/src/basis/basdeclsP.h new file mode 100644 index 0000000..7527fe3 --- /dev/null +++ b/src/basis/basdeclsP.h @@ -0,0 +1,39 @@ +#ifndef _BASDECLSP_H +#define _BASDECLSP_H + +// declarations for substitution by cpp for compressed array count meanings +// only for capitalized versions + +// define HEAD_NTAGS 1 +// define HEAD_NCONT 2 +// define HEAD_NPRIM 3 +// define HEAD_NCOEF 4 +// define TAG_NCONT 1 +// define TAG_NPRIM 2 +// define TAG_NCOEF 3 +// define TAG_FCONT 4 +// define TAG_LCONT 5 +// define CONT_TYPE 1 +// define CONT_NPRIM 2 +// define CONT_NGEN 3 +// define CONT_IEXP 4 +// define CONT_ICFP 5 +// define CONT_TAG 6 + +#define HEAD_NTAGS 1 +#define HEAD_NCONT 2 +#define HEAD_NPRIM 3 +#define HEAD_NCOEF 4 +#define TAG_NCONT 1 +#define TAG_NPRIM 2 +#define TAG_NCOEF 3 +#define TAG_FCONT 4 +#define TAG_LCONT 5 +#define CONT_TYPE 1 +#define CONT_NPRIM 2 +#define CONT_NGEN 3 +#define CONT_IEXP 4 +#define CONT_ICFP 5 +#define CONT_TAG 6 + +#endif \ No newline at end of file diff --git a/src/basis/basis.c b/src/basis/basis.c new file mode 100644 index 0000000..e2fc025 --- /dev/null +++ b/src/basis/basis.c @@ -0,0 +1,1410 @@ +#include +#include +#include + +#include "inp.h" +#include "basP.h" +#include "geom.h" +#include "rtdb.h" +#include "geomP.h" +#include "context.h" +#include "basdeclsP.h" +#include "geobasmapP.h" + +/* + Block data structure to initialize the common block variables in the + internal basis set object data structures +*/ +int nbasis_rtdb = 0; +bool bsactive[nbasis_bsmx] = {false}; // Assuming nbasis_bsmx is a predefined constant +int angular_bs[nbasis_bsmx] = {-565}; // Assuming nbasis_bsmx is a predefined constant +double bsversion = 1.00; + +bool bas_version() { +/* + Routine that calclulates the size of the common block structures + used in the basis set object and the mapped representation object. + input none + output always true. +*/ + int cdata, idata, rdata; + int mapidata, total4, total8; + + // character data + cdata = 256 * 2 * nbasis_bsmx + 256 * nbasis_rtdb_mx; + cdata = cdata + 16 * ntags_bsmx * nbasis_bsmx; + + // real data + rdata = mxbs_exndcf * nbasis_bsmx + 1; + rdata = 8 * rdata; + + // integer data in basis set object common + idata = ndbs_head * nbasis_bsmx; + idata = idata + ndbs_tags * ntags_bsmx * nbasis_bsmx; + idata = idata + ndbs_ucont * nucont_bsmx * nbasis_bsmx; + idata = idata + nbasis_bsmx * 4 + 1; + idata = idata + nbasis_rtdb_mx; + idata = 4 * idata; + + // integer data in the mapped object. + mapidata = 4 * ncont_mx * nbasis_bsmx; + mapidata = mapidata + 3 * nat_mx * nbasis_bsmx; + mapidata = mapidata + 4 * nbasis_bsmx; + mapidata = 4 * mapidata; + + // total space + total4 = idata + mapidata; + total8 = 2 * total4 + rdata + cdata; + total4 = total4 + rdata + cdata; + + printf(" **** basis set version %f ****\n", bsversion); + printf(" character data in-core %d bytes\n", cdata); + printf(" real data in-core %d bytes\n", rdata); + printf(" integer*4 data in-core %d bytes\n", idata); + printf("or integer*8 data in-core %d bytes\n", (2 * idata)); + printf(" integer*4 mapping data in-core %d bytes\n", mapidata); + printf("or integer*8 mapping data in-core %d bytes\n", (2 * mapidata)); + printf(" total(4) = %d bytes\n", total4); + printf(" total(8) = %d bytes\n", total8); + + // convert to kilobytes + cdata = (cdata + 999) / 1000; + rdata = (rdata + 999) / 1000; + idata = (idata + 999) / 1000; + mapidata = (mapidata + 999) / 1000; + total4 = (total4 + 999) / 1000; + total8 = (total8 + 999) / 1000; + + printf(" **** basis set version %f ****\n", bsversion); + printf(" character data in-core %d Kbytes\n", cdata); + printf(" real data in-core %d Kbytes\n", rdata); + printf(" integer*4 data in-core %d Kbytes\n", idata); + printf("or integer*8 data in-core %d Kbytes\n", (2 * idata)); + printf(" integer*4 mapping data in-core %d Kbytes\n", mapidata); + printf("or integer*8 mapping data in-core %d Kbytes\n", (2 * mapidata)); + printf(" total(4) = %d Kbytes\n", total4); + printf(" total(8) = %d Kbytes\n", total8); + + return true; +} + +bool bas_create(FILE *basis, char *name) { +/* + creates a handle and marks it active in the in-core data structure + + passed + integer basis [output] returned handle + char*(*) name [input] name of basis set. +*/ + int i; + + for (i = 0; i <= nbasis_bsmx; i++) { + if (!bsactive[i]) break; + } + + if (i > nbasis_bsmx) { + printf(" bas_create: no free basis handles for %s\n", name); + return false; + } + + // store some information in basis data structure + // (NOTE: name discarded in LOAD operation) + bs_name[basis] = name; + len_bs_name[basis] = strlen(name); + + // Initialize basis info to be empty + bs_trans[basis] = " "; + int *infbs_head[basis] = malloc(ndbs_head * sizeof(int)); + double *exndcf[basis] = malloc(mxbs_exndcf * sizeof(double)); + int *infbs_tags[basis] = malloc(ndbs_tags * ntags_bsmx * sizeof(int)); + int *infbs_cont[basis] = malloc(ndbs_ucont * nucont_bsmx * sizeof(int)); + + // Initialize geo-basis info to empty + int *ibs_cn2ucn[basis] = malloc( ncont_mx * sizeof(int)); + int *ibs_cn2ce[basis] = malloc(ncont_mx * sizeof(int)); + int *ibs_ce2uce[basis] = malloc(nat_mx * sizeof(int)); + int *ibs_cn2bfr[basis] = malloc(2 * ncont_mx * sizeof(int)); + int *ibs_ce2cnr[basis] = malloc(2 * nat_mx * sizeof(int)); + ncont_tot_gb[basis] = 0; + nprim_tot_gb[basis] = 0; + nbf_tot_gb[basis] = 0; + ibs_geom[basis] = 0; + + // Mark basis as active and return info + bsactive[basis] = true; + basis = basis - BASIS_HANDLE_OFFSET; + return true; + +} + +bool bas_destroy(FILE *basisin) { + /* + destroys information about an active incore basis + and the associated mapping arrays. + */ + + bool ret_val; + int basis; + + ret_val = bas_check_handle(basisin, "bas_destroy"); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + + if (!gbs_map_clear(basisin)) { + printf(" error clearing map "); + return false; + } + + bsactive[basis] = false; + + return true; +} + +bool bas_check_handle(FILE *basisin, char *msg) { + + // Checks to see if a basis handle is valid + + // local variables + bool ret_val; + int basis; + basis = basisin + BASIS_HANDLE_OFFSET; + ret_val = (basis > 0 && basis <= nbasis_bsmx); + if (ret_val) { + ret_val = ret_val && bsactive[basis]; + } + + if (!ret_val) { + printf("%s: basis handle is invalid\n", msg); + printf("basis_check_handle: lexical handle %d\n", basis); + printf("basis_check_handle: handle %d\n", basisin); + } + return ret_val; +} + +bool bas_321g_load(FILE *rtdb) { + /* + routine to load the rtdb with 321g basis sets for atoms 1-20 + H to Ca + */ + + // local variables + bool ret_val; + int basis, usebas, itag, icont; + + if (!bas_create(basis, "321g:1-20")) { + printf("basis set handle not created\n"); + return false; + } + + usebas = basis + BASIS_HANDLE_OFFSET; + bs_name[usebas] = "321g:1-20"; + + // only h and o for now + + // hydrogen and oxygen + bs_tags[0][usebas] = 'H'; + bs_tags[1][usebas] = 'He'; + bs_tags[2][usebas] = 'O'; + infbs_head[HEAD_NTAGS][usebas] = 3; + infbs_head[HEAD_NCONT][usebas] = 9; + infbs_head[HEAD_NPRIM][usebas] = 30; + infbs_head[HEAD_NCOEF][usebas] = 30; + + // hydrogen + itag = 0; + infbs_tags[TAG_NCONT][itag][usebas] = 2; + infbs_tags[TAG_NPRIM][itag][usebas] = 3; + infbs_tags[TAG_NCOEF][itag][usebas] = 3; + infbs_tags[TAG_FCONT][itag][usebas] = 1; + infbs_tags[TAG_LCONT][itag][usebas] = 2; + + // hydrogen cont 1 bas 1 + icont = 0; + infbs_cont[CONT_TYPE][icont][usebas] = 0; + infbs_cont[CONT_NPRIM][icont][usebas] = 2; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 1; + infbs_cont[CONT_ICFP][icont][usebas] = 3; + infbs_cont[CONT_TAG][icont][usebas] = itag; + exndcf[1][usebas] = 5.44717800; + exndcf[2][usebas] = 0.82454700; + exndcf[3][usebas] = 0.15628500; + exndcf[4][usebas] = 0.90469100; + + // hydrogen cont 2 bas 2 + icont++; + infbs_cont[CONT_TYPE][icont][usebas] = 0; + infbs_cont[CONT_NPRIM][icont][usebas] = 1; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 5; + infbs_cont[CONT_ICFP][icont][usebas] = 6; + infbs_cont[CONT_TAG][icont][usebas] = itag; + exndcf[5][usebas] = 0.18319200; + exndcf[6][usebas] = 1.00000000; + + // helium + itag++; + icont++; + infbs_tags[TAG_NCONT][itag][usebas] = 2; + infbs_tags[TAG_NPRIM][itag][usebas] = 3; + infbs_tags[TAG_NCOEF][itag][usebas] = 3; + infbs_tags[TAG_FCONT][itag][usebas] = icont; + infbs_tags[TAG_LCONT][itag][usebas] = icont + 2 - 1; + + // helium cont 1 bas 3 + infbs_cont[CONT_TYPE][icont][usebas] = 0; + infbs_cont[CONT_NPRIM][icont][usebas] = 2; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 7; + infbs_cont[CONT_ICFP][icont][usebas] = 9; + infbs_cont[CONT_TAG][icont][usebas] = itag; + exndcf[7][usebas] = 13.62670000; + exndcf[8][usebas] = 1.99935000; + exndcf[9][usebas] = 0.17523000; + exndcf[10][usebas] = 0.89348300; + + // helium cont 2 bas 4 + icont++; + infbs_cont[CONT_TYPE][icont][usebas] = 0; + infbs_cont[CONT_NPRIM][icont][usebas] = 1; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 11; + infbs_cont[CONT_ICFP][icont][usebas] = 12; + infbs_cont[CONT_TAG][icont][usebas] = itag; + exndcf[11][usebas] = 0.38299300; + exndcf[12][usebas] = 1.00000000; + + // oxygen + itag++; + icont++; + infbs_tags[TAG_NCONT][itag][usebas] = 5; + infbs_tags[TAG_NPRIM][itag][usebas] = 9; + infbs_tags[TAG_NCOEF][itag][usebas] = 9; + infbs_tags[TAG_FCONT][itag][usebas] = icont; + infbs_tags[TAG_LCONT][itag][usebas] = icont + 5 - 1; + + // oxygen cont 1 bas 5 + infbs_cont[CONT_TYPE][icont][usebas] = 0; + infbs_cont[CONT_NPRIM][icont][usebas] = 3; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 13; + infbs_cont[CONT_ICFP][icont][usebas] = 16; + infbs_cont[CONT_TAG][icont][usebas] = itag; + exndcf[13][usebas] = 322.03700000; + exndcf[14][usebas] = 48.43080000; + exndcf[15][usebas] = 10.42060000; + exndcf[16][usebas] = 0.05923940; + exndcf[17][usebas] = 0.35150000; + exndcf[18][usebas] = 0.70765800; + + // oxygen cont 2 bas 6 + icont++; + infbs_cont[CONT_TYPE][icont][usebas] = 0; + infbs_cont[CONT_NPRIM][icont][usebas] = 2; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 19; + infbs_cont[CONT_ICFP][icont][usebas] = 21; + infbs_cont[CONT_TAG][icont][usebas] = 2; + exndcf[19][usebas] = 7.40294000; + exndcf[20][usebas] = 1.57620000; + exndcf[21][usebas] = -0.40445300; + exndcf[22][usebas] = 1.22156000; + + // oxygen cont 3 bas 7 + icont++; + infbs_cont[CONT_TYPE][icont][usebas] = 1; + infbs_cont[CONT_NPRIM][icont][usebas] = 2; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 23; + infbs_cont[CONT_ICFP][icont][usebas] = 25; + infbs_cont[CONT_TAG][icont][usebas] = 2; + exndcf[23][usebas] = 7.40294000; + exndcf[24][usebas] = 1.57620000; + exndcf[25][usebas] = 0.24458600; + exndcf[26][usebas] = 0.85395500; + + // oxygen cont 4 bas 8 + icont++; + infbs_cont[CONT_TYPE][icont][usebas] = 0; + infbs_cont[CONT_NPRIM][icont][usebas] = 1; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 27; + infbs_cont[CONT_ICFP][icont][usebas] = 28; + infbs_cont[CONT_TAG][icont][usebas] = 2; + exndcf[27][usebas] = 0.37368400; + exndcf[28][usebas] = 1.00000000; + + // oxygen cont 5 bas 9 + icont++; + infbs_cont[CONT_TYPE][icont][usebas] = 1; + infbs_cont[CONT_NPRIM][icont][usebas] = 1; + infbs_cont[CONT_NGEN][icont][usebas] = 1; + infbs_cont[CONT_IEXP][icont][usebas] = 29; + infbs_cont[CONT_ICFP][icont][usebas] = 30; + infbs_cont[CONT_TAG][icont][usebas] = 2; + exndcf[29][usebas] = 0.37368400; + exndcf[30][usebas] = 1.00000000; + + ret_val = bas_rtdb_do_store(rtdb, bs_name[usebas], bs_tags[1][usebas], + infbs_head[1][usebas], infbs_tags[1][1][usebas], + infbs_cont[1][1][usebas], exndcf[1][usebas], 3, 9, 30); + + printf("inside 321g load\n"); + printf("inside 321g load\n"); + ret_val = bas_print(basis); + printf("inside 321g load\n"); + printf("inside 321g load\n"); + + if (!bas_destroy(basis)) { + printf("error releasing temporary basis handle\n"); + return false; + } + + return true; +} + +bool bas_print(FILE *basisin) { + // routine to print unique basis information that is in core + + int mytags, myucont, myprim, mycoef, basis, len_tag, empty; + int i, j, k, l, ifcont, mygen, mytype, iexptr, icfptr; + char ctype[7][4] = {"S", "P", "D", "F", "G", "H", "I"}; + char cltype[3][4] = {"SP", "SPD"}; + char shell_type[4]; + char blank[17] = " "; + bool ret_val; + + basis = basisin + BASIS_HANDLE_OFFSET; + + ret_val = bas_check_handle(basisin, "bas_print"); + if (!ret_val) return false; + + // print basis set information + printf(" Basis \"%.*s\" -> \"%.*s\"\n", inp_strlen(bs_name[basis]), bs_name[basis], inp_strlen(bs_trans[basis]), bs_trans[basis]); + + mytags = infbs_head[HEAD_NTAGS][basis]; + if (mytags <= 0) { + printf(" Empty basis set \n\n"); + return true; + } + + myucont = infbs_head[HEAD_NCONT][basis]; + myprim = infbs_head[HEAD_NPRIM][basis]; + mycoef = infbs_head[HEAD_NCOEF][basis]; + + /* write(6,2) mytags, myucont, myprim, mycoef + 2 format( + $ ' number of unique tags :',i5/ + $ ' number of unique contractions :',i5/ + $ ' number of unique primitives :',i5/ + $ ' number of unique coefficients :',i5) + */ + + for (i = 1; i <= mytags; i++) { + len_tag = inp_strlen(bs_tags[i][basis]); + empty = (16 - len_tag) / 2; + printf("%.*s%s\n", empty, blank, bs_tags[i][basis]); + + myucont = infbs_tags[TAG_NCONT][i][basis]; +/* write(6,4) myucont,infbs_tags(TAG_NPRIM,i,basis), + $ infbs_tags(TAG_NCOEF,i,basis) + 4 format( + & ' number of contractions:',i5/ + & ' number of primitives :',i5/ + & ' number of coefficients:',i5/) +*/ + ifcont = infbs_tags[TAG_FCONT][i][basis]; + + printf(" Exponent Coefficients \n"); + printf(" ----------- %.*s\n", 60, "----------------"); + + for (j = 1; j <= myucont; j++) { + myprim = infbs_cont[CONT_NPRIM][ifcont][basis]; + mygen = infbs_cont[CONT_NGEN][ifcont][basis]; + + mytype = infbs_cont[CONT_TYPE][ifcont][basis]; + if (mytype < 0) { + strcpy(shell_type, cltype[abs(mytype)]); + } else { + strcpy(shell_type, ctype[mytype]); + } + +/* write(6,5) j, shell_type(1:inp_strlen(shell_type)), +* $ myprim, mygen +* 5 format(/ +* $ ' contraction :', i5/ +* $ ' type : ',a/ +* $ ' number of primitives :', i5/ +* $ ' number of contractions:', i5/) +*/ + iexptr = infbs_cont[CONT_IEXP][ifcont][basis] - 1; + icfptr = infbs_cont[CONT_ICFP][ifcont][basis] - 1; + + for (k = 1; k <= myprim; k++) { + printf(" %2d %2.*s %11.4f", j, 2, shell_type, exndcf[iexptr + k][basis]); + for (l = 1; l <= mygen; l++) { + printf(" %10.6f", exndcf[icfptr + k + (l - 1) * myprim][basis]); + } + printf("\n"); + } + printf("\n"); + + ifcont++; + } + } + +/* If geom is set print out the info about total basis info associated with + the geometry also + + ... not done yet +*/ + + return true; +} + +bool bas_rtdb_load(FILE *rtdb, FILE *geom, FILE *basisin, char *name) { +/* + arguments + rtdb -> valid rtdb handle + geom -> valid geometry handle with info loaded + basisin -> valid basis handle + name -> basis set name that must exist on the rtdb + passed +*/ + int basis; + bool ret_val; + int lentmp; + char tmp[256]; + bool rtdb_status, status; + int i, nat, idum_cont, idum_at; + int j, jstart, jend, jsize; + int kstart, kend, ksize, lsize, icount; + int nbf, iu_cont, myang; + bool foundit; + + rtdb_status = true; + + // check geom and basis handles returns false if either is invalid + + ret_val = geom_check_handle(geom, "bas_rtdb_load"); + if (!ret_val) return false; + ret_val = bas_check_handle(basisin, "bas_rtdb_load"); + if (!ret_val) return false; + + // store geom tag with basis map info + basis = basisin + BASIS_HANDLE_OFFSET; + ibs_geom[basis] = geom; + + // translate "name" to current "context" + bs_name[basis] = name; + len_bs_name[basis] = strlen(name); + if (!context_rtdb_match(rtdb, name, bs_trans[basis])) + bs_trans[basis] = name; + len_bs_trans[basis] = strlen(bs_trans[basis]); + + // generate rtdb names and load information + strcpy(tmp, "basis:"); + strncat(tmp, bs_trans[basis], len_bs_trans[basis]); + lentmp = strlen(tmp) + 1; + + strcpy(tmp + lentmp, ":bs_tags"); + rtdb_status = rtdb_status && rtdb_par_cget(rtdb, tmp, ntags_bsmx, bs_tags[1][basis]); + + strcpy(tmp + lentmp, ":exps and coeffs"); + rtdb_status = rtdb_status && rtdb_par_get(rtdb, tmp, mt_dbl, mxbs_exndcf, exndcf[1][basis]); + + strcpy(tmp + lentmp, ":header"); + rtdb_status = rtdb_status && rtdb_par_get(rtdb, tmp, mt_int, ndbs_head, infbs_head[1][basis]); + + strcpy(tmp + lentmp, ":tags info"); + rtdb_status = rtdb_status && rtdb_par_get(rtdb, tmp, mt_int, ndbs_tags * ntags_bsmx, infbs_tags[1][1][basis]); + + strcpy(tmp + lentmp, ":contraction info"); + rtdb_status = rtdb_status && rtdb_par_get(rtdb, tmp, mt_int, ndbs_ucont * nucont_bsmx, infbs_cont[1][1][basis]); + + // read the basis now get check status of read operations + if (!rtdb_status) { + printf("bas_rtdb_load: ERROR\n"); + printf("basis set is not there\n"); + printf("name requested <|%s|>\n", bs_name[basis]); + printf("translated name requested <|%s|>\n", bs_trans[basis]); + return false; + } + + status = geom_ncent(geom, &nat); + if (nat == 0 || !status) { + printf("bas_rtdb_load: ERROR\n"); + printf("number of centers is zero or weird\n"); + printf("nat = %d\n", nat); + return false; + } +/* + build center to unique center map + + do 00100 i = 1,nat + if (.not.inp_match(infbs_head(HEAD_NTAGS,basis),.true., + & tags(i,geom), + & bs_tags(1,basis),ibs_ce2uce(i,basis))) then + write(6,*)' geom tag was not found in basis tag list rtdb' + write(6,*)' geom tag searched for ',tags(i,geom) + write(6,*)' basis tag list was from basis ',bs_name(basis) + c.... add further diagnostics later + endif + #if defined(BS_DEBUG) + write(6,*)' ibs_ce2uce(',i,') = ',ibs_ce2uce(i,basis) + #endif + 00100 continue +*/ + // build center to unique center map + for (i = 0; i <= nat; i++) { + foundit = false; + for (j = 0; j <= infbs_head[HEAD_NTAGS][basis]; j++) { + if (inp_compare(true, tags[i][geom], bs_tags[j][basis])) { + ibs_ce2uce[i][basis] = j; + foundit = true; + break; + } + } + if (!foundit) { + printf("geom tag was not found in basis tag list rtdb\n"); + printf("geom tag searched for %s\n", tags[i][geom]); + printf("basis tag list was from basis %s\n", bs_name[basis]); + } + } + + // build total # of contractions and center -> contraction range map + idum_cont = infbs_tags[TAG_NCONT][ibs_ce2uce[0][basis]][basis]; + ncont_tot_gb[basis] = idum_cont; + ibs_ce2cnr[0][0][basis] = 1; + ibs_ce2cnr[1][0][basis] = idum_cont; + for (i = 1; i <= nat; i++) { + idum_cont = infbs_tags[TAG_NCONT][ibs_ce2uce[i][basis]][basis]; + ncont_tot_gb[basis] = idum_cont + ncont_tot_gb[basis]; + ibs_ce2cnr[0][i][basis] = ibs_ce2cnr[1][i - 1][basis] + 1; + ibs_ce2cnr[1][i][basis] = ibs_ce2cnr[0][i][basis] + idum_cont - 1; + } + + // build contraction -> center map + for (i = 0; i <= nat; i++) { + jstart = ibs_ce2cnr[0][i][basis]; + jend = ibs_ce2cnr[1][i][basis]; + for (j = jstart; j <= jend; j++) { + ibs_cn2ce[j][basis] = i; + } + } + + // build contraction -> unique contraction map + for (i = 0; i <= nat; i++) { + idum_at = ibs_ce2uce[i][basis]; + jstart = ibs_ce2cnr[0][i][basis]; + jend = ibs_ce2cnr[1][i][basis]; + jsize = jend - jstart + 1; + kstart = infbs_tags[TAG_FCONT][idum_at][basis]; + kend = infbs_tags[TAG_LCONT][idum_at][basis]; + ksize = kend - kstart + 1; + lsize = infbs_tags[TAG_NCONT][idum_at][basis]; + if (jsize == ksize && ksize == lsize) { + icount = 0; + for (j = jstart; j <= jend; j++) { + ibs_cn2ucn[j][basis] = kstart + icount; + icount++; + } + } else { + printf("bas_rtdb_load: ERROR\n"); + printf("contraction range size <-> unique contraction range size mismatch\n"); + printf(" contraction range (%d:%d)\n", jstart, jend); + printf(" unique contraction range (%d:%d)\n", kstart, kend); + printf(" contraction size: %d\n", jsize); + printf(" calculated unique contraction size: %d\n", ksize); + printf(" lookup unique contraction size: %d\n", lsize); + return false; + } + } + + // build nprim_tot_gb, nbf_tot_gb, and + // contraction -> basis function range map + + iu_cont = ibs_cn2ucn[0][basis]; + nbf = infbs_cont[CONT_TYPE][iu_cont][basis]; + nbf = (nbf + 1) * (nbf + 2) / 2; + nbf = nbf * infbs_cont[CONT_NGEN][iu_cont][basis]; + + ibs_cn2bfr[0][0][basis] = 1; + ibs_cn2bfr[1][0][basis] = nbf; + + nbf_tot_gb[basis] = nbf; + nprim_tot_gb[basis] = infbs_cont[CONT_NPRIM][iu_cont][basis]; + for (i = 1; i <= ncont_tot_gb[basis]; i++) { + iu_cont = ibs_cn2ucn[i][basis]; + + nbf = infbs_cont[CONT_TYPE][iu_cont][basis]; + nbf = (nbf + 1) * (nbf + 2) / 2; + nbf = nbf * infbs_cont[CONT_NGEN][iu_cont][basis]; + + ibs_cn2bfr[0][i][basis] = ibs_cn2bfr[1][i - 1][basis] + 1; + ibs_cn2bfr[1][i][basis] = ibs_cn2bfr[0][i][basis] + nbf - 1; + + nbf_tot_gb[basis] += nbf; + nprim_tot_gb[basis] += infbs_cont[CONT_NPRIM][iu_cont][basis]; + } + + // build high angular momentum of this loaded pair + + for (i = 0; i <= ncont_tot_gb[basis]; i++) { + iu_cont = ibs_cn2ucn[i][basis]; + myang = infbs_cont[CONT_TYPE][iu_cont][basis]; + angular_bs[basis] = max(angular_bs[basis], myang); + } + + return ret_val; +} + +bool bas_rtdb_store(FILE *rtdb, char *name, FILE *basisin) { + + /* + Store basis set (not geometry) related info about specified + basis in into the rtdb with the given name + */ + + bool ret_val; + int basis; + + int ret_val = bas_check_handle(basisin, "bas_rtdb_store"); + if (!ret_val) return false; + basis = basisin + BASIS_HANDLE_OFFSET; + ret_val = bas_rtdb_do_store(rtdb, name, + bs_tags[0][basis], infbs_head[0][basis], + infbs_tags[0][0][basis], + infbs_cont[0][0][basis], exndcf[basis], + infbs_head[HEAD_NTAGS][basis], + infbs_head[HEAD_NCONT][basis], + infbs_head[HEAD_NPRIM][basis] + infbs_head[HEAD_NCOEF][basis]); + + return ret_val; +} + +bool bas_rtdb_do_store(FILE *rtdb, char *name, int *tagsin, int *head_array, + int *tags_array, int *ucont_array, double *excfin, + int ntagsin, int nucontin, int nexcf) { +/* + This routine stores the basis set information in the appropriate + data structure on the run-time-data-base (rtdb). + + This is a private routine called by the user level routine + bas_rtdb_store(rtdb, name, basis) + + argument description + C*(*) name -> name of the basis set + C*16(1:ntagsin) tagsin -> name of each tag +*/ + bool status; + char tmp[256]; + int len_name, lentmp; + + bool bas_rtdb_do_store = true; + + status = bas_rtdb_in(rtdb); + + // generate rtdb names and store information + len_name = strlen(name); + strcpy(tmp, "basis:"); + strncat(tmp, name, len_name); + lentmp = strlen(tmp) + 1; + + status = true; + status = status && bas_rtdb_add(rtdb, name); + strncpy(tmp + lentmp, ":bs_tags", 9); + status = status && rtdb_par_cput(rtdb, tmp, ntagsin, tagsin); + + strncpy(tmp + lentmp, ":exps and coeffs", 16); + status = status && rtdb_par_put(rtdb, tmp, mt_dbl, nexcf, excfin); + + strncpy(tmp + lentmp, ":header", 8); + status = status && rtdb_par_put(rtdb, tmp, mt_int, ndbs_head, head_array); + + strncpy(tmp + lentmp, ":tags info", 10); + status = status && rtdb_par_put(rtdb, tmp, mt_int, (ndbs_tags * ntagsin), tags_array); + + strncpy(tmp + lentmp, ":contraction info", 16); + status = status && rtdb_par_put(rtdb, tmp, mt_int, (ndbs_ucont * nucontin), ucont_array); + + // read the basis now and check status of read operations + if (!status) { + printf("bas_rtdb_store: ERROR\n"); + printf("one or more put operations failed\n"); + bas_rtdb_do_store = false; + // add diagnostics later + return bas_rtdb_do_store; + } + + return ret_val; +} + +bool bas_high_angular(FILE *basisin, int *high_angular) { + +// calculate and store high angular momentem function +// for given basis. + + bool ret_val; + int basis, myucont, i; + + ret_val = bas_check_handle(basisin, "bas_high_angular"); + if (!ret_val) { + printf("basis handle not valid\n"); + return false; + } + + basis = basisin + BASIS_HANDLE_OFFSET; + if (angular_bs[basis] > -565) { + *high_angular = angular_bs[basis]; + return true; + } + + myucont = infbs_head[HEAD_NCONT][basis]; + *high_angular = -565; + + for (i = 1; i <= myucont; i++) { + *high_angular = max(*high_angular, infbs_cont[CONT_TYPE][i][basis]); + } + + angular_bs[basis] = high_angular; + return true; +} + +bool gbs_map_clear(FILE *basisin) { + +// routine to clear online map information and basis information + + int basis; + + basis = basisin + BASIS_HANDLE_OFFSET; + + if (!bas_check_handle(basisin, "gbs_map_clear")) { + printf("basis handle not valid\n"); + return false; + } + + ifill(ncont_mx, 0, ibs_cn2ucn(1, basis), 1); + ifill(ncont_mx, 0, ibs_cn2ce(1, basis), 1); + ifill(2 * ncont_mx, 0, ibs_cn2bfr(1, 1, basis), 1); + ifill(nat_mx, 0, ibs_ce2uce(1, basis), 1); + ifill(2 * nat_mx, 0, ibs_ce2cnr(1, 1, basis), 1); + ncont_tot_gb[basis] = 0; + nprim_tot_gb[basis] = 0; + nbf_tot_gb[basis] = 0; + + angular_bs[basis] = -565; + + return true; +} + +bool gbs_map_print(FILE *basisin) { + + FILE *mygeom; + int nat, basis, i, myfirst, mylast, mysize, mycenter, myucont; + bool status, ret_val; + + basis = basisin + BASIS_HANDLE_OFFSET; + + // check geom and basis handles returns false if either is invalid + mygeom = ibs_geom(basis); + ret_val = geom_check_handle(mygeom, "gbs_map_print"); + if (!ret_val) return false; + ret_val = bas_check_handle(basisin, "gbs_map_print"); + if (!ret_val) return false; + + // find number of atoms + status = geom_ncent(mygeom, &nat); + + if (nat == 0 || !status) { + printf("gbs_map_print: ERROR\n"); + printf("number of centers is zero or weird\n"); + printf("nat = %d\n", nat); + ret_val = false; + // add diagnostics later + return ret_val; + } + + // print global information + printf("<<< GBS_MAP_PRINT >>>\n"); + printf("total number of atoms : %d\n", nat); + printf("total number of contractions : %d\n", ncont_tot_gb[basis]); + printf("total number of primitives : %d\n", nprim_tot_gb[basis]); + printf("total number of basis functions : %d\n", nbf_tot_gb[basis]); + + // print center based mapping information + printf("\n"); + printf("=============================================================\n"); + printf("center -> unique center map \n"); + printf(" -> contraction range map \n"); + printf("=============================================================\n"); + for (i = 1; i <= nat; i++) { + printf("center: %d maps to unique center: %d\n", i, ibs_ce2uce(i, basis)); + myfirst = ibs_ce2cnr(1, i, basis); + mylast = ibs_ce2cnr(2, i, basis); + mysize = mylast - myfirst + 1; + printf("has %d contractions \n", mysize, myfirst, mylast); + } + + // print contraction based mapping information + printf("\n"); + printf("=============================================================\n"); + printf("contraction -> center map \n"); + printf(" -> unique contraction in basis set \n"); + printf(" -> basis function range \n"); + printf("=============================================================\n"); + + for (i = 1; i <= ncont_tot_gb[basis]; i++) { + mycenter = ibs_cn2ce(i, basis); + myucont = ibs_cn2ucn(i, basis); + myfirst = ibs_cn2bfr(1, i, basis); + mylast = ibs_cn2bfr(2, i, basis); + mysize = mylast - myfirst + 1; + printf("contraction %d is on center: %d\n", i, mycenter); + printf("is represented by unique contraction: %d\n", myucont); + printf("has %d basis functions \n", mysize, myfirst, mylast); + } + + return true; +} + +bool bas_get_exponent(FILE *basisin, int icont, bool unique, double *exp) { + // passed + bool ret_val; + int basis, myucont, icontmax; + int myprim, myexptr; + + ret_val = bas_check_handle(basisin, "bas_get_exponent"); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + if (unique) { + icontmax = infbs_head(HEAD_NCONT, basis); + myucont = icont; + } else { + icontmax = ncont_tot_gb(basis); + myucont = ibs_cn2ucn(icont, basis); + } + ret_val = icont > 0 && icont <= icontmax; + if (!ret_val) { + printf("bas_get_exponent: ERROR\n"); + printf("contraction range for basis is 1:%d\n", icontmax); + printf("information requested for contraction:%d\n", icont); + return false; + } + + myexptr = infbs_cont(CONT_IEXP, myucont, basis); + myprim = infbs_cont(CONT_NPRIM, myucont, basis); + dcopy(myprim, exndcf(myexptr, basis), 1, exp, 1); + /* + for (int i = 0; i < myprim; i++) { + exp[i] = exndcf(myexptr, basis)[i]; + } + */ + + return true; +} + +bool bas_continfo(FILE *basisin, int icont, bool unique, int *nprimo, int *ngeno, int *sphcart) { + + // passed + bool ret_val; + int basis, myucont, icontmax; + + nprimo = -123; + ngeno = -456; + sphcart = -789; + + ret_val = bas_check_handle(basisin, "bas_continfo"); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + + if (unique) { + icontmax = infbs_head(HEAD_NCONT, basis); + } else { + icontmax = ncont_tot_gb(basis); + } + + if (!(icont > 0 && icont <= icontmax)) { + printf("bas_continfo: ERROR\n"); + if (unique) { + printf("unique contraction range for basis is 1:%d\n", icontmax); + } else { + printf("contraction range for basis is 1:%d\n", icontmax); + } + printf("information requested for contraction:%d\n", icont); + return false; + } + + if (unique) { + myucont = icont; + } else { + myucont = ibs_cn2ucn(icont, basis); + } + + // ... no spherical yet 3/94 only cart. + sphcart = 0; + nprimo = infbs_cont(CONT_NPRIM, myucont, basis); + ngeno = infbs_cont(CONT_NGEN, myucont, basis); + + return true; +} + +bool bas_numcont(FILE *basisin, int *numcont, bool unique) { + int basis; + numcont = -6589; + if (!bas_check_handle(basisin, "bas_numcont")) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + + if (unique) { + numcont = infbs_head(HEAD_NCONT, basis); + } else { + numcont = ncont_tot_gb(basis); + } + + return true; +} + +bool bas_numbf(FILE *basisin, int *nbf) { + int basis; + *nbf = -6589; + if (!bas_check_handle(basisin, "bas_numbf")) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + *nbf = nbf_tot_gb(basis); + return true; +} + +bool bas_get_coeff(FILE *basisin, int icont, bool unique, double *coeff) { + bool ret_val; + int basis, myucont, icontmax; + int mycoeffptr, myprim, mygen; + + ret_val = bas_check_handle(basisin, "bas_get_coeff"); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + + if (unique) { + icontmax = infbs_head(HEAD_NCONT, basis); + myucont = icont; + } else { + icontmax = ncont_tot_gb(basis); + myucont = ibs_cn2ucn(icont, basis); + } + + ret_val = (icont > 0) && (icont <= icontmax); + if (!ret_val) { + printf("bas_get_coeff: ERROR\n"); + printf("contraction range for basis is 1:%d\n", icontmax); + printf("information requested for contraction:%d\n", icont); + return false; + } + + mycoeffptr = infbs_cont(CONT_ICFP, myucont, basis); + myprim = infbs_cont(CONT_NPRIM, myucont, basis); + mygen = infbs_cont(CONT_NGEN, myucont, basis); + dcopy((myprim * mygen), exndcf(mycoeffptr, basis), 1, coeff, 1); + + return true; +} + +bool bas_set_exponent(FILE *basisin, int icont, bool unique, double *exp, int nexp) { + bool ret_val; + int basis, myucont, icontmax; + int myexptr, myprim; + + ret_val = bas_check_handle(basisin, "bas_set_exponent"); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + + if (unique) { + icontmax = infbs_head(HEAD_NCONT, basis); + myucont = icont; + } else { + icontmax = ncont_tot_gb(basis); + myucont = ibs_cn2ucn(icont, basis); + } + + ret_val = (icont > 0) && (icont <= icontmax); + if (!ret_val) { + printf("bas_set_exponent: ERROR\n"); + printf("contraction range for basis is 1:%d\n", icontmax); + printf("information requested for contraction:%d\n", icont); + return false; + } + + myexptr = infbs_cont(CONT_IEXP, myucont, basis); + myprim = infbs_cont(CONT_NPRIM, myucont, basis); + ret_val = (myprim == nexp); + if (!ret_val) { + printf("bas_set_exponent: ERROR\n"); + printf("input and stored number of exponents (nprim) differ\n"); + printf("input nprim: %d\n", nexp); + printf("stored nprim: %d\n", myprim); + return false; + } + + dcopy(nexp, exp, 1, exndcf(myexptr, basis), 1); + + return true; +} + +bool bas_set_coeff(FILE *basisin, int icont, bool unique, double *coeff, int ncoeff) { + bool ret_val; + int basis, myucont, icontmax; + int mycoeffptr, myprim, mygen; + + ret_val = bas_check_handle(basisin, "bas_set_coeff"); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + + if (unique) { + icontmax = infbs_head(HEAD_NCONT, basis); + myucont = icont; + } else { + icontmax = ncont_tot_gb(basis); + myucont = ibs_cn2ucn(icont, basis); + } + + ret_val = (icont > 0) && (icont <= icontmax); + if (!ret_val) { + printf("bas_set_coeff: ERROR\n"); + printf("contraction range for basis is 1:%d\n", icontmax); + printf("information requested for contraction:%d\n", icont); + return false; + } + + mycoeffptr = infbs_cont(CONT_ICFP, myucont, basis); + myprim = infbs_cont(CONT_NPRIM, myucont, basis); + mygen = infbs_cont(CONT_NGEN, myucont, basis); + + ret_val = (ncoeff == (myprim * mygen)); + if (!ret_val) { + printf("bas_set_coeff: ERROR\n"); + printf("input and stored number of coefficients (nprim*ngen) differ\n"); + printf("input nprim*ngen: %d\n", ncoeff); + printf("stored nprim*ngen: %d\n", (myprim * mygen)); + return false; + } + + dcopy(ncoeff, coeff, 1, exndcf(mycoeffptr, basis), 1); + + return true; +} + +bool bas_rtdb_out(FILE *rtdb) { + bool ret_val; + + ret_val = rtdb_par_put(rtdb, "basis:nbasis", MT_INT, 1, nbasis_rtdb) && + rtdb_par_cput(rtdb, "basis:names", nbasis_rtdb, bs_names_rtdb); + if (!ret_val) { + printf(" bas_rtdb_out: rtdb is corrupt "); + } + + return bas_rtdb_out; +} + +bool bas_rtdb_add(FILE *rtdb, char *name) { + int basis; + bool ret_val, status; + int ln; + + // See if name is on the rtdb already + ln = strlen(name); + status = bas_rtdb_in(rtdb); + ret_val = true; + + for (basis = 1; basis <= nbasis_rtdb; basis++) { + if (strncmp(name, bs_names_rtdb[basis], ln) == 0) { + return ret_val; + } + } + + // Name is not present ... add and rewrite info + if (nbasis_rtdb == nbasis_rtdb_mx) { + printf("bas_rtdb_add: too many basetries on rtdb %s\n", name); + return false; + } + + nbasis_rtdb++; + strcpy(bs_names_rtdb[nbasis_rtdb], name); + len_bs_rtdb[nbasis_rtdb] = ln; + + ret_val = bas_rtdb_out(rtdb); + if (!ret_val) { + printf("bas_rtdb_add: rtdb error adding %.*s\n", ln, name); + return ret_val; + } + + return true; +} + +bool bas_print_all() { + bool ret_val; + int basis, basin; + + ret_val = true; + for (basis = 1; basis <= nbasis_bsmx; basis++) { + if (bsactive(basis)) { + basin = basis - BASIS_HANDLE_OFFSET; + ret_val = ret_val && bas_print(basin); + } + } + + return ret_val; +} + +void bas_err_info(char *info) { + int bas, basin; + int nbas; + bool status; + + /* + For internal use of the basis set routines only: print out + info of known basis sets to aid in diagnosing a problem + */ + nbas = 0; + for (bas = 1; bas <= nbasis_bsmx; bas++) { + if (bsactive(bas)) { + nbas++; + } + } + printf("%s: open basis sets: %d\n", info, nbas); + + nbas = 0; + for (bas = 1; bas <= nbasis_bsmx; bas++) { + if (bsactive(bas)) { + basin = bas - BASIS_HANDLE_OFFSET; + status = bas_print(basin); + } + } + + if (nbasis_rtdb > 0) { + printf("%s: basis sets in current rtdb %d\n", info, nbasis_rtdb); + for (bas = 1; bas <= nbasis_rtdb; bas++) { + printf("number: %d basis set name: %.*s\n", bas, len_bs_rtdb[bas], bs_names_rtdb[bas]); + } + } +} + +bool bas_rtdb_in(FILE *rtdb) { + bool ret_val; + int bas; + /* + load in info about known basis sets ... this is more + for diagnostic and debugging purposes + */ + ret_val = false; + nbasis_rtdb = 0; + if (rtdb_par_get(rtdb, "basis:nbasis", MT_INT, 1, &nbasis_rtdb)) { + if (!rtdb_par_cget(rtdb, "basis:names", nbasis_rtdb_mx, bs_names_rtdb)) { + printf("bas_rtdb_in: rtdb corrupt\n"); + } else { + for (bas = 1; bas <= nbasis_rtdb; bas++) { + len_bs_rtdb[bas] = inp_strlen(bs_names_rtdb[bas]); + } + ret_val = true; + } + } + + return ret_val; +} + +bool bas_cn2ce(FILE *basisin, int cont, int *center) { + bool ret_val; + int basis; + + ret_val = bas_check_handle(basisin); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + ret_val = cont > 0 && cont <= ncont_tot_gb(basis); + if (!ret_val) { + printf("bas_cn2ce: invalid contraction information\n"); + printf("contraction range is 1:%d\n", ncont_tot_gb(basis)); + printf("input contraction was: %d\n", cont); + return false; + } + center = ibs_cn2ce(cont, basis); + + return true; +} + +bool bas_cn2bf(FILE *basisin, int cont, int *ifirst, int *ilast) { + bool ret_val; + int basis; + + ret_val = bas_check_handle(basisin); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + ret_val = cont > 0 && cont <= ncont_tot_gb(basis); + if (!ret_val) { + printf("bas_cn2bf: invalid contraction information\n"); + printf("contraction range is 1:%d\n", ncont_tot_gb(basis)); + printf("input contraction was: %d\n", cont); + return false; + } + + ifirst = ibs_cn2bfr(1, cont, basis); + ilast = ibs_cn2bfr(2, cont, basis); + + return true; +} + +bool bas_ce2cnr(FILE *basisin, int center, int ifirst, int ilast) { + bool ret_val; + int basis, nat; + + ret_val = bas_check_handle(basisin); + if (!ret_val) return false; + + basis = basisin + BASIS_HANDLE_OFFSET; + ret_val = geom_ncent(ibs_geom(basis), nat); + if (nat == 0 || !ret_val) { + printf("bas_ce2cnr: ERROR\n"); + printf("number of centers is zero or weird\n"); + printf("nat = %d\n", nat); + // add diagnostics later + return false; + } + + ret_val = center > 0 && center <= nat; + if (!ret_val) { + printf("bas_ce2cnr: invalid center information\n"); + printf("contraction range is 1:%d\n", nat); + printf("input contraction was: %d\n", center); + return false; + } + + return true; +} + +bool bas_add_ucnt(FILE *basis, char *tag, int l_value, int ngen, int nprim, + double *expnt, double *coeffs, int ldc) { + bool ret_val; + int ind; // Index into basis function structures + int free; // Free space pointer + int i, itag, jtag, iu_cont, ntags; // Locals + + /* + adds a new general contraction on the specified tag. If the + tag is not present it will also add that by calling bas_add_utag + */ + + ret_val = bas_check_handle(basis, "bas_add_ucnt"); + if (!ret_val) return false; + ind = basis + BASIS_HANDLE_OFFSET; + + // Make sure that the tag is in the list + ret_val = bas_add_utag(basis, tag, &itag); + if (!ret_val) return false; + + // Update header information about all unique contractions on all tags + free = infbs_head[HEAD_NPRIM][ind] + infbs_head[HEAD_NCOEF][ind] + 1; + if ((free+nprim*ngen+nprim-1) > mxbs_exndcf) { + printf("bas_add_ucnt: too many prims/coeffs\n"); + return false; + } + if (infbs_head[HEAD_NCONT][ind]+1 > nucont_bsmx) { + printf("bas_add_ucnt: too many contractions\n"); + return false; + } + + infbs_head[HEAD_NCONT][ind] = infbs_head[HEAD_NCONT][ind] + 1; + infbs_head[HEAD_NPRIM][ind] = infbs_head[HEAD_NPRIM][ind] + nprim; + infbs_head[HEAD_NCOEF][ind] = infbs_head[HEAD_NCOEF][ind] + ngen*nprim; + + ntags = infbs_head[HEAD_NTAGS][ind]; + if (itag != ntags) { + for (jtag = ntags; jtag >= itag+1; jtag--) { + // Shuffle data+pointers for following tags up one contraction + for (iu_cont = infbs_tags[TAG_LCONT][jtag][ind]; + iu_cont >= infbs_tags[TAG_FCONT][jtag][ind]; iu_cont--) { + for (i = 1; i <= ndbs_ucont; i++) { + infbs_cont[i][iu_cont+1][ind] = infbs_cont[i][iu_cont][ind]; + } + } + // Increment first and last contractions on following tags + infbs_tags[TAG_FCONT][jtag][ind] = infbs_tags[TAG_FCONT][jtag][ind] + 1; + infbs_tags[TAG_LCONT][jtag][ind] = infbs_tags[TAG_LCONT][jtag][ind] + 1; + } + } + + // Increment basis info on this tag + infbs_tags[TAG_NCONT][itag][ind] = infbs_tags[TAG_NCONT][itag][ind] + 1; + infbs_tags[TAG_NPRIM][itag][ind] = infbs_tags[TAG_NPRIM][itag][ind] + nprim; + infbs_tags[TAG_NCOEF][itag][ind] = infbs_tags[TAG_NCOEF][itag][ind] + nprim*ngen; + if (infbs_tags[TAG_FCONT][itag][ind] == 0) { + if (itag != ntags) { + printf("bas_add_ucnt: tag error %d\n", itag); + } + infbs_tags[TAG_FCONT][itag][ind] = infbs_head[HEAD_NCONT][ind]; + infbs_tags[TAG_LCONT][itag][ind] = infbs_head[HEAD_NCONT][ind]; + } else { + infbs_tags[TAG_LCONT][itag][ind] = infbs_tags[TAG_LCONT][itag][ind] + 1; + } + + iu_cont = infbs_tags[TAG_LCONT][itag][ind]; // Index of new contraction + + infbs_cont[CONT_TYPE][iu_cont][ind] = l_value; + infbs_cont[CONT_NPRIM][iu_cont][ind] = nprim; + infbs_cont[CONT_NGEN][iu_cont][ind] = ngen; + infbs_cont[CONT_TYPE][iu_cont][ind] = l_value; + infbs_cont[CONT_IEXP][iu_cont][ind] = free; + infbs_cont[CONT_ICFP][iu_cont][ind] = free + nprim; + + // Copy real data over + memcpy(&exndcf[free][ind], expnt, nprim * sizeof(double)); + free = free + nprim; + for (i = 1; i <= ngen; i++) { + memcpy(&exndcf[free][ind], &coeffs[(i-1)*ldc], nprim * sizeof(double)); + free = free + nprim; + } + /* + write(6,*) ' expnt input ', nprim + call output(expnt, 1, nprim, 1, 1, nprim, 1, 1) + write(6,*) ' coeffs input ', nprim, ngen, ldc + call output(coeffs, 1, nprim, 1, ngen, ldc, ngen, 1) + + Done + */ + + return true; +} + +bool bas_add_utag(FILE *basisin, char *tag, int *itag) { + bool ret_val; + int basis; // [local] index into basis arrays + + // Add the unique tag to the list of tags in the basis, + // incrementing the no. of tags if necessary. + // Return in itag the index of the unique tag + + ret_val = bas_check_handle(basisin, "bas_add_utag"); + if (!ret_val) return false; + basis = basisin + BASIS_HANDLE_OFFSET; + + for (*itag = 1; *itag <= infbs_head[HEAD_NTAGS][basis]; (*itag)++) { + if (inp_compare(true, bs_tags[*itag][basis], tag)) return true; + } + + // No match found ... append new tag to the list + + *itag = infbs_head[HEAD_NTAGS][basis] + 1; + if (*itag > ntags_bsmx) { + printf("bas_add_utag: too many tags %d\n", itag); + return false; + } + + infbs_head[HEAD_NTAGS][basis] = itag; + bs_tags[itag][basis] = tag; + + return true; +} diff --git a/src/basis/doc/api b/src/basis/doc/api new file mode 100644 index 0000000..a6d6b59 --- /dev/null +++ b/src/basis/doc/api @@ -0,0 +1,138 @@ + + + + Program main + + ---------------------- + + fix rtdb for arrays of characters + + + ---------------------- + + context management + + + ---------------------- + + get/set coords,charges,tags,masses,zmat(?),ncenters, + + map 'geometry' -> name of geometry + + logical geom_load(rtdb, 'geometry', geom) + logical geom_store(rtdb, 'geometry', geom) + + ncent = geom_ncenter(geom) + call geom_tag(geom, icent, tag) + call geom_cent_coords(geom, icent, coords) + call geom_cent_charge(geom, icent, charge) + ... + call geom_cent_info(geom, icent, tag, coords, charge, mass, ...) + + logical geom_zmat_defined() + + nvariables = geom_zmat_nvars(geom) + nconstants .... + + call geom_cart_get(geom, all info) + call geom_cart_set(geom, all info) + + call geom_zmat_get + call geom_zmat_set + + print + + + Also on the DB + + - list of known geometry names + + + ---------------------- + + map 'mo basis' -> name of basis descriptor + + ---------------------- + + logical basis_load(name_of_basis_descriptor, geom, basis) + nbasis_func = basis_nfunc(basis) + nbasis_shell = basis_nshell(basis) + natoms / basis_centers = + map atom/center<->shell<->bf + get/set exponents/contraction coeffs + shell info (angular, gcontract, spherical/cart) + highest ang. mom. + print + load/store + + On the data-base is + + - list of known basis set names + - + + ---------------------- + + Cannot tweak geometry or basis between init/term calls + +... control + int_initialize(geom, num_basis, basis_array) : generate internal int structures +................................................................................. +# not needed if batmol writes to rtdb properly + int_initialize_tape10 () : generate internal int structures (batmol?) +................................................................................. + int_terminate() : throw away internal int structures + int_print_known_basis() + int_set_eri_timing() + int_report_eri_timings() + int_mem(max1e, maxg, mscratch_1e, mscratch_2e) + int_mem_one(max1e, mscratch_1e) + int_mem_4(maxg, mscratch_2e) + int_mem_3(max3, mscratch_3_2e) + int_mem_2(max2, mscratch_2_2e) + int_mem_3ov(max3ov, mscratch_3ov) + +... two electron +. 4 center 2e integrals + eri = + int_two_4(bra_basis, ket_basis, ish, jsh, ksh, lsh, lscr, + scr, eri) + lab_two_4(bra_basis, ket_basis, ish, jsh, ksh, lsh, zerotol, + canonicalize, eri, nints, ilab, jlab, klab, llab) + +. 3 center 2e integrals + eri = + int_two_3 (bra_basis, ket_basis, ish, jsh, ksh, lscr, scr, eri) + lab_two_3 (bra_basis, ket_basis, canonical_bra, canonical_both, ish, jsh, ksh, zerotol, + eri, nints, ilab, jlab, klab) + +. 2 center 2e integrals + eri = + int_two_2 (bra_basis, ket_basis, ish, jsh, lscr, scr, eri) + lab_two_2 (bra_basis, ket_basis, canonical_both, ish, jsh, zerotol, lscr, scr, eri, + nints, ilab, jlab) + +... one electron integrals + int_one_ke_basic (i_basis, j_basis, ish, jsh, lscr, scr, T) + int_one_pe_basic (i_basis, j_basis, ish, jsh, lscr, scr, V) + int_one_ov_basic (i_basis, j_basis, ish, jsh, lscr, scr, S) + int_one_h1_basic (i_basis, j_basis, ish, jsh, lscr, scr, H1) + int_one_all_basic(i_basis, j_basis, ish, jsh, lscr, scr, S, T, V) + lab_one_ke (i_basis, j_basis, ish, jsh, zerotol, ilab, jlab, T, numt) + lab_one_pe (i_basis, j_basis, ish, jsh, zerotol, ilab, jlab, V, numv) + lab_one_ov (i_basis, j_basis, ish, jsh, zerotol, ilab, jlab, S, nums) + lab_one_h1 (i_basis, j_basis, ish, jsh, zerotol, ilab, jlab, H1, numh1) + lab_one_all(i_basis, j_basis, ish, jsh, zerotol, ilab, jlab, S, T, V, numstv) + + one_3c_int = + int_one_3ov(i_basis, j_basis, k_basis, ish, jsh, ksh, lscr, + scr, OV3) + lab_one_3ov(i_basis, j_basis, k_basis, ish, jsh, ksh, zerotol, + OV3, ilab, jlab, klab, numov3) + + + int_mpole(i_basis, j_basis, Lvalue, ish, jsh, lscr, scr, MPINTS) + lab_mpole(i_basis, j_basis, Lvalue, ish, jsh, MPINTS, ilab, jlab, zerotol) + + + periodic versions (with k vector) + ---------------------- + \ No newline at end of file diff --git a/src/basis/doc/basis.doc b/src/basis/doc/basis.doc new file mode 100644 index 0000000..e9fef9f --- /dev/null +++ b/src/basis/doc/basis.doc @@ -0,0 +1,149 @@ +/* The basis set objects are written in C. + +Proposal for basis set objects (rjh/rak) + + + 1) Currently have only general and segmented contractions of + primitive gaussians but other basis sets should be anticipated + + 2) The whole GTO basis is either cartesian or spherical harmonic + + 3) Basis functions are associated with atomic tags, not coordinates, + the tags providing the connection to a geometry + + 4) All basis functions associated with an 'atomic' center will be + numbered consecutively + +5) General Basis Set Class: + +**** Attributes: +Basis_type: Contracted Gaussian + or pseudo-potentails + or plane wave + or .?.?.?.?. + +if (Basis_type.eq.1) (sub class definition) + number_of_tags: number of tags with information to be given/retrieved. + (NOTE: generally equal to the number of centers to be + defined by the geometry object.) + cartesian: is it cartesian or spherical (transformed). + nprim_tot:: total number of primitive gaussians. + ncoeff_tot:: total number of contractraction coeffs + nshell_tot:: total number of shells (in the normal sense). + nbf_tot:: total number of basis functions (in the normal sense). + + + + for each (tag (1 .. number_of_tags) (a.k.a. center)) + tag: name of atomic center (o, bond_function, ghost center) + contraction_type: Segmented or General? + (note: basis sets may have mixed segmented and + general contracted centers) + ngen: =1 (segmented) >1 (general). + nshell: number of shells for this tag (in normal sense). + :nGshell: number of generally contracted shells. + (what is the limit?? tradeoff memory vs cpu) + nbf_on_tag: number of basis functions for this tag. + :nbf_cart: number of cartesian basis functions + (note: evaluation most likely always in cartesians) + :nbf_sph: number of spherical basis functions. + nprim: number of primitives in contraction. + + for each (contraction set) [general or segmented] + type basis function type = GTO (redundant??) + L angular momentum s,p,d = 0, 1, 2 + Ltype =0 use only angular momentum L + =1 implies use of all angular momenta up to L + e.g., sp shells (like gaussian) or + spd shells (Rydberg, bond functions) + num_cnt number of contractions coefs (always=1 for segmented) + nbf number of basis functions in contraction (external view) + :nbf_cart: number of cartesian basis functions in contraction + :nbf_sph: number of spherical basis functions in contraction + nprim number of primitives in contraction set + coef(nprim,num_cnt) contraction coefs + ex(nprim) exponents + + +proposal: (some of what we need to decide!!) +. above atributes that are :attribute: are to be used only by the api + and integral routines. +. attributes that are atribute:: are derived some how. + + +other pointer arrays etc can be defined/derived from this information +and the goemetry object and thus are application dependent. There is +no requirement that these pointer arrays match across applications/ +modules/libraries. Mappings to the final integral code could be +supplied by the api. The major hash will be to determine what is +defined a priori and what is derived information. e.g., the total +number of basis functions for a basis: is fixed but for a +molecular or periodic system it is unknown until the geometry object +defines the scope. + +else + it ain't been defined yet. +endif + +**** Operations: (lean and mean) +query_on: scope_set, number_of_tags, cartesian, nprim_tot, ncoeff_tot, nbf_tot etc. +open: +close: +define: (later not now). +set_scope: + + + 6) For compact and portable representation in the database this is + compacted into five entries. + + a) dimension information + + b) character information + + c) pointer information + + d) integer information + + e) real information + + 7) All basis set info is named as if components of a basis module + so all names are of the form + + basis:: ... + + 8) A data base entry basis:list contains a list of the names + of all basis sets in the database so that it is easy to + examine what is in the database + + 9) a mixed basis set can be a sum of two basis sets object subclasses?? + e.g., planewave + contracted gaussian + pseudo-potential + contracted gaussian. + + + +objects: + basis set --> name --> on rtdb +operations: + ---------------------- + + map 'mo basis' -> name of basis descriptor + + ---------------------- + + logical basis_load(name_of_basis_descriptor, geom, basis) + nbasis_func = basis_nfunc(basis) + nbasis_shell = basis_nshell(basis) + natoms / basis_centers = + map atom/center<->shell<->bf + get/set exponents/contraction coeffs + shell info (angular, gcontract, spherical/cart) + highest ang. mom. + print + load/store + + On the data-base is + + - list of known basis set names + - + +*/ \ No newline at end of file diff --git a/src/basis/doc/basis.output b/src/basis/doc/basis.output new file mode 100755 index 0000000..9c9b9ac --- /dev/null +++ b/src/basis/doc/basis.output @@ -0,0 +1,139 @@ + SI 0 + S 20 1.00 + 3948000.00000000 0.00000204 + 591100.00000000 0.00001584 + 134500.00000000 0.00008336 + 38120.00000000 0.00035136 + 12460.00000000 0.00127660 + 4504.00000000 0.00415191 + 1758.00000000 0.01230300 + 729.10000000 0.03331020 + 318.00000000 0.08098450 + 144.60000000 0.17029000 + 67.97000000 0.28687900 + 32.82000000 0.33034000 + 16.03000000 0.19660200 + 7.39600000 0.03545350 + 3.66100000 -0.00053520 + 1.82300000 0.00161465 + 0.91470000 -0.00037274 + 0.33930000 0.00014623 + 0.15000000 -0.00007894 + 0.06438000 0.00001928 + S 20 1.00 + 3948000.00000000 -0.00000054 + 591100.00000000 -0.00000422 + 134500.00000000 -0.00002218 + 38120.00000000 -0.00009360 + 12460.00000000 -0.00034012 + 4504.00000000 -0.00111061 + 1758.00000000 -0.00330878 + 729.10000000 -0.00911602 + 318.00000000 -0.02287900 + 144.60000000 -0.05171190 + 67.97000000 -0.09990910 + 32.82000000 -0.15274700 + 16.03000000 -0.12750800 + 7.39600000 0.09469630 + 3.66100000 0.41403600 + 1.82300000 0.46793400 + 0.91470000 0.17392700 + 0.33930000 0.00843895 + 0.15000000 -0.00099807 + 0.06438000 0.00036210 + S 20 1.00 + 3948000.00000000 0.00000014 + 591100.00000000 0.00000108 + 134500.00000000 0.00000569 + 38120.00000000 0.00002395 + 12460.00000000 0.00008724 + 4504.00000000 0.00028416 + 1758.00000000 0.00084984 + 729.10000000 0.00233527 + 318.00000000 0.00590466 + 144.60000000 0.01334610 + 67.97000000 0.02628890 + 32.82000000 0.04074260 + 16.03000000 0.03614760 + 7.39600000 -0.03039230 + 3.66100000 -0.13596100 + 1.82300000 -0.25014400 + 0.91470000 -0.15805000 + 0.33930000 0.36965500 + 0.15000000 0.61771800 + 0.06438000 0.22251400 + S 1 1.00 + 0.91470000 1.00000000 + S 1 1.00 + 0.33930000 1.00000000 + S 1 1.00 + 0.15000000 1.00000000 + S 1 1.00 + 0.06438000 1.00000000 + P 12 1.00 + 1780.00000000 0.00020121 + 421.80000000 0.00174937 + 136.70000000 0.00948141 + 51.81000000 0.03723130 + 21.60000000 0.11076300 + 9.56300000 0.23793300 + 4.35000000 0.35369100 + 2.00600000 0.32883900 + 0.92050000 0.13237300 + 0.35000000 0.01033000 + 0.13810000 -0.00015031 + 0.05338000 0.00026581 + P 12 1.00 + 1780.00000000 -0.00004272 + 421.80000000 -0.00037704 + 136.70000000 -0.00202240 + 51.81000000 -0.00812833 + 21.60000000 -0.02422720 + 9.56300000 -0.05438250 + 4.35000000 -0.07990510 + 2.00600000 -0.08889580 + 0.92050000 0.01839970 + 0.35000000 0.33509600 + 0.13810000 0.53228800 + 0.05338000 0.25437400 + P 1 1.00 + 0.92050000 1.00000000 + P 1 1.00 + 0.35000000 1.00000000 + P 1 1.00 + 0.13810000 1.00000000 + P 1 1.00 + 0.05338000 1.00000000 + D 1 1.00 + 0.12600000 1.00000000 + D 1 1.00 + 0.32100000 1.00000000 + D 1 1.00 + 0.81700000 1.00000000 + D 1 1.00 + 2.08200000 1.00000000 + F 1 1.00 + 0.16900000 1.00000000 + F 1 1.00 + 0.34100000 1.00000000 + F 1 1.00 + 0.68800000 1.00000000 + G 1 1.00 + 0.32000000 1.00000000 + G 1 1.00 + 0.70500000 1.00000000 + H 1 1.00 + 0.58300000 1.00000000 + S 1 1.00 + 0.02600000 1.00000000 + P 1 1.00 + 0.01920000 1.00000000 + D 1 1.00 + 0.04680000 1.00000000 + F 1 1.00 + 0.07350000 1.00000000 + G 1 1.00 + 0.15100000 1.00000000 + H 1 1.00 + 0.32300000 1.00000000 + **** \ No newline at end of file diff --git a/src/basis/doc/robert.doc b/src/basis/doc/robert.doc new file mode 100644 index 0000000..783d62e --- /dev/null +++ b/src/basis/doc/robert.doc @@ -0,0 +1,218 @@ + + +1) Minimize implementation effort + +2) Simplify data structures so that are flattened + more readily + +3) Enable local integral routines to work directly + from API interface and/or internal data structures + +4) Store info in the database to avoid having a + zillion small files floating around + + +In the database we store just the basis set description +(i.e., the atomic basis sets for the unique atom tags) +in as simple a format as possible + +In core, we have in addition mapping arrays that build +the basis set from the geometry and the basis set description + + +How to store the atomic basis set compactly, but so that +it is readily stored and efficiently used? + + The integral routines will be given basis set handles and + shell indices. These will be used to lookup + the shell info (l-value, ngen, nprim) and find pointers + to the contraction info. Since Fortran cannot have pointers + returned to it from a C interface we have to either store + the stuff both on the C and Fortran sides, or do it + all in Fortran (I know, double ugh). I see no point + in doing things twice. + + Since we are stuck with F77 we have no structures and + are back using simple offsets etc. This actually makes + storing the info externally easier since the internal + representation is flat. + + +Detailed data structures ... derive from their usage by your +integral routines. + + int_2e_4c(ibasis, jbasis, ish, jsh, ksh, lsh, ...) + + check basis handles + + get info (type, nprim, ngen, coords) on each shell + + find pointers to coeffs/exponents for each shell + (this implies that they are stored packed into a + single array and we have offsets stored) + + branch to the fastest routine depending on if generally contracted, + the angular momentum, if it is an sp shell, ... + + in your API + + call the primitive evaluation routine with explicit + coord/coeff/exponents + + + So it seems that we have very similar data structures to the + present int.h, except that the basis info is only stored for + unique atom types + + Now do the mapping in detail + + + if (ibasis .le.0 .or. ibasis .gt. nbasis) call errquit(...) + + if (ish .le. 0 .or. ish .gt. nshell(ibasis)) call errquit(...) + + iuniq = shell_uniq(ish, ibasis) ... map shell to no. of the shell info + for unique tags only + + itype = shell_type(iuniq, ibasis) (1, 2, 3 for s, p, d + -1, -2, ... for sp, spd, ... shells) + + iprim = shell_nprim(iuniq, ibasis) + + igen = shell_nprim(iuniq, ibasis) + + iexpnt= shell_expt(iuniq, ibasis) ... offset in exp(1, ibasis) where + this shells exponents start + + icoeff= shell_cofpt(iuniq, ibasis) ... offset in coeff(1, ibasis) .... + + + icent= shell_cent(ish, ibasis) ... center no. for this shell to get coords + (for efficiency should grab the + coords from the geometry) + + + Should be rolling at this point. + + Also need the following arrays to support the other basis set + routines + + cent_to_sh(1:2, icent, ibasis) (contains hi-lo) + cent_to_bf(1:2, icent, ibasis) + sh_to_bf(1:2, ish, ibasis) + + + The info about each shell is simply + + integer type, nprim, ngen + real coeff(nprim,ngen), expnt(nprim) + + + Thus, the atomic basis set is just + + integer nshell, nprim_tag, ncoeff_tag + integer type(nshell), nprim(nshell), ngen(nshell), + cofpt(nshell), expt(nsehll) + real coeff(ncoeff_tag), expnt(nprim_tag) + + + For external storage this can be compactly represented as + ... and there is no reason why the data cannot be also used + this way (so that coefpt and expt provide offsets into rdata) + + integer dim_info(3) + integer idata(5*nshell) + real rdata(ncoeff_tag+nprim_tag) + + (the rtdb can automatically allocate the MA arrays and read into them) + + These could be stored on the rtdb as + + basis:basis_name:tag:dim_info + basis:basis_name:tag:idata + basis:basis_name:tag:rdata + + along with a summary of all unique tag info + + integer nshell_total, nprim_total, idata_total, rdata_total + + basis:basis_name:dim_info -> integer dim_info(4) + + + + However, we can make things even easier by storing the whole + damn lot in one data structure since it will always be possible + to store info on the unqiue atom centers (even if the whole + periodic table is in there!). Thus, my recomendation is that + the data base contain the following + + basis:basis_name:dim_info integer + basis:basis_name:tags character + basis:basis_name:tdata integer + basis:basis_name:idata integer + basis:basis_name:rdata double precision + + Where + + dim_info(1) -> nshell_uniq_total = total no. of shells on the unique tags + dim_info(2) -> nprim_uniq_total = total no. of prims on the unique tags + dim_info(3) -> idata_uniq_total = total length of idata + dim_info(4) -> idata_uniq_total = total length of rdata + dim_info(5) -> ntags_uniq = no. of unique tags + + tags(1:ntags_uniq) = character array of tags (cannot be allocated + using MA !!) + + tdata(1, itag_uniq) = first unique shell on this tag + tdata(2, itag_uniq) = last last shell on this tag + + idata(1, ish_uniq) = type of shell + idata(2, ish_uniq) = nprim in shell + idata(3, ish_uniq) = ngen of shell + idata(4, ish_uniq) = offset into rdata for coeffs + idata(5, ish_uniq) = offset into rdata for exponents + idata(6, ish_uniq) = no. of bf in this shell + + + + To load this lot into core and build the data structures + on the fly : + + logical function basis_load(rtdb, name, igeom, ibasis) + + 0) look for translations of name within the current or higher + context using context_rtdb_match(). With the name or available + translation look for basis:basis_name:dim_info ... if this is there + then the basis set is defined. Can adopt a default at this + point if desired. Check that have statically allocated enuf + space to read in the unique tags. + + rtdb_cget( tags ) + rtdb_ma_get (tdata, rdata, idata) + + 1) Get tags/coords info from the geometry (note .. only one geometry + being used by the integrals at a time ... I would suggest that + the geometry handle be removed from the int_init() call and + be stored internal to each basis sets structure ... this then + gives us a mechanism to compute integrals between different + geometries (this sounds worth thinking about more)). + + 2) Loop thru centers checking that have a basis defined for that + tag and accumulate the no. of shells and basis functions. + At same time build map from atoms to shells and bf and + map from shells to unique shell no. + + +Done. + +How this info gets onto the database is another problem. I would +suggest that the input program provides the info in nearly this form +to a basis set routine for output to the database. Since the input +routines want to be very general it's best to let them worry about +the details and live with a very simple basis set interface. We also +need to consider how to handle plane waves, giaos etc ... let's talk +about this and also to Jeff about giaos before casting this in +FORTRAN. I think that the above will suffice for the GTO basis +sets and we can add additional RTDB entries for the plane waves etc. + +Robert \ No newline at end of file diff --git a/src/basis/geobasmapP.h b/src/basis/geobasmapP.h new file mode 100644 index 0000000..b6c5b97 --- /dev/null +++ b/src/basis/geobasmapP.h @@ -0,0 +1,29 @@ +#ifndef _GEOBASMAPP_H +#define _GEOBASMAPP_H +/* + NOTE: this MUST follow basP.fh in the include order + + these are all in core + mapping arrays atoms <-> contr <-> bfn (maybe add shells later) + + contraction -> unique_contr :: ibs_cn2ucn(ncont,nbasis) + contraction -> center_number :: ibs_cn2ce (ncont,nbasis) + center -> unique_center :: ibs_ce2uce(nat,nbasis) + contraction -> basis function range :: ibs_cn2bfr(2,ncont,nbasis) + centers -> contraction range :: ibs_ce2cnr(2,nat,nbasis) +*/ +#ifdef __cplusplus +extern "C" { +#endif + int ibs_cn2ucn(int ncont, int nbasis); + int ibs_cn2ce(int ncont, int nbasis); + int ibs_ce2uce(int nat, int nbasis); + int ibs_cn2bfr(int ncont, int nbasis); + int ibs_ce2cnr(int nat, int nbasis); + + int ncont_tot_gb, nprim_tot_gb, nbf_tot_gb, ibs_geom; +#ifdef __cplusplus +} +#endif + +#endif // _GEOBASMAPP_H diff --git a/src/basis/newbasis.c b/src/basis/newbasis.c new file mode 100644 index 0000000..96869f5 --- /dev/null +++ b/src/basis/newbasis.c @@ -0,0 +1,11 @@ +#include + +void bas_dummmmmm() { + int i; + i = 0; +} + +int bas_print_known(FILE *rtdb) { + return 0; +} + diff --git a/src/basis/testbasis.c b/src/basis/testbasis.c new file mode 100644 index 0000000..d9190ab --- /dev/null +++ b/src/basis/testbasis.c @@ -0,0 +1,170 @@ +#include "bas.h" +#include "rtdb.h" +#include "geom.h" + +#include +#include + +int main() { + FILE *rtdb, *geom, *basis; + int ngen, nprim, iang; + int ncenters, sphcart, i, j; + char drivtags[20][16]; + double coords[3][20], charge[20]; + double exp[400], coeff[400]; + bool status; + double expnt_new[3] = {1.0, 2.0, 3.0}; + double coeff_new[4][3] = { + {-1.0, -2.0, -3.0}, + {0.0, -4.0, -5.0}, + {-6.0, 0.0, -7.0}, + {-8.0, -9.0, 0.0} + }; + + if (!ma_init(MT_DBL, -1, -1)) { + printf("Error initializing ma_init\n"); + return 99; + } + + status = rtdb_par_open("shit.rtdb", "unknown", &rtdb); + + printf("rtdb handle %d\n", rtdb); + + status = bas_321g_load(&rtdb); + + if (!geom_create(&geom, "321g:1-20")) { + printf("Error getting geometry handle\n"); + return 1; + } + + ncenters = 3; + // oxygen + strcpy(drivtags[0], "O"); + coords[0][0] = 0.0; + coords[1][0] = 0.0; + coords[2][0] = 0.0; + charge[0] = 8.0; + // hydrogen 1 + strcpy(drivtags[1], "H"); + coords[0][1] = 1.0; + coords[1][1] = 1.0; + coords[2][1] = 1.0; + charge[1] = 1.0; + // hydrogen 2 + strcpy(drivtags[2], "H"); + coords[0][2] = -1.0; + coords[1][2] = -1.0; + coords[2][2] = -1.0; + charge[2] = 1.0; + + if (!geom_cart_set(&geom, ncenters, drivtags, coords, charge)) { + printf("geom_cart_set fail\n"); + return 1; + } else { + status = geom_print(&geom); + } + + basis = 0; + if (!bas_create(&basis, "3211-20")) { + printf("Error getting basis handle\n"); + return 1; + } + + basis = 0; + if (!bas_create(&basis, "321g1-20")) { + printf("Error getting second basis handle\n"); + return 1; + } + + basis = 0; + if (!bas_create(&basis, "321g:1-20")) { + printf("Error getting third basis handle\n"); + return 1; + } + + status = bas_rtdb_load(&rtdb, &geom, &basis, "321g:1-20"); + status = bas_print(&basis); + status = gbs_map_print(&basis); + + status = bas_continfo(&basis, 1, false, &nprim, &ngen, &sphcart); + printf("f:query: nprim cont 1 %d\n", nprim); + printf("f:query: ngen cont 1 %d\n", ngen); + printf("f:query: sphcart cont 1 %d\n", sphcart); + + status = bas_get_exponent(&basis, 1, false, exp); + status = bas_get_coeff(&basis, 1, false, coeff); + printf("exponents and coefficients\n"); + for (i = 0; i < nprim; i++) { + for (j = 0; j < ngen; j++) { + printf("%lf ", coeff[i + j * nprim]); + } + printf("\n"); + } + + status = bas_continfo(&basis, 1, true, &nprim, &ngen, &sphcart); + printf("t:query: nprim cont 1 %d\n", nprim); + printf("t:query: ngen cont 1 %d\n", ngen); + printf("t:query: sphcart cont 1 %d\n", sphcart); + + status = bas_get_exponent(&basis, 1, true, exp); + status = bas_get_coeff(&basis, 1, true, coeff); + printf("exponents and coefficients\n"); + for (i = 0; i < nprim; i++) { + for (j = 0; j < ngen; j++) { + printf("%lf ", coeff[i + j * nprim]); + } + printf("\n"); + } + + exp[0] = 565.6589; + coeff[0] = 6.021023; + status = bas_set_exponent(&basis, 1, true, exp, nprim + 1); + status = bas_set_coeff(&basis, 1, true, coeff, nprim * ngen + 1); + status = bas_set_exponent(&basis, 1, true, exp, nprim); + status = bas_set_coeff(&basis, 1, true, coeff, nprim * ngen); + + status = bas_continfo(&basis, 1, true, &nprim, &ngen, &sphcart); + printf("modified\n"); + printf("t:query: nprim cont 1 %d\n", nprim); + printf("t:query: ngen cont 1 %d\n", ngen); + printf("t:query: sphcart cont 1 %d\n", sphcart); + + status = bas_get_exponent(&basis, 1, true, exp); + status = bas_get_coeff(&basis, 1, true, coeff); + printf("exponents and coefficients\n"); + for (i = 0; i < nprim; i++) { + for (j = 0; j < ngen; j++) { + printf("%lf ", coeff[i + j * nprim]); + } + printf("\n"); + } + + // Try adding new contractions on an existing center + printf("adding 3*3 d function on H\n"); + if (!bas_add_ucnt(basis, "H", 2, 3, 3, expnt_new, coeff_new, 4)) { + printf(" basis_add_ucnt failed"); + } + if (!bas_print(&basis)) { + printf(" print ???"); + } + + // Try adding new contractions on a new center + printf("adding 2*3 g function on Cl\n"); + if (!bas_add_ucnt(basis, "Cl", 4, 2, 3, expnt_new, coeff_new, 4)) { + printf(" basis_add_ucnt failed"); + } + if (!bas_print(&basis)) { + printf(" print ???"); + } + + printf("bas_print_all\n"); + status = bas_print_all(); + bas_err_info("who who who"); + + status = bas_high_angular(&basis, &iang); + printf("high angular momentum %d\n", iang); + status = bas_version(); + + printf("testbasis done\n"); + return 0; +} diff --git a/src/config/makefile.h b/src/config/makefile.h new file mode 100644 index 0000000..1669144 --- /dev/null +++ b/src/config/makefile.h @@ -0,0 +1,176 @@ + +# $Id: makefile.h,v 1.1.1.1 1994-03-29 06:44:24 d3g681 Exp $ + +# Common definitions for all makefiles ... these can be overridden +# either in each makefile by putting additional definitions below the +# include statement, or on the command line + + +# +# Set TOPDIR to point to your top-level directory that contains +# src, lib, config, ... (SRCDIR, etc., are derived from TOPDIR) +# + TOPDIR = /msrc/home/d3g681/allnew + SRCDIR = $(TOPDIR)/src + LIBDIR = $(TOPDIR)/lib + BINDIR = $(TOPDIR)/bin + INCDIR = $(TOPDIR)/src/include + +# +# Define TARGET to be the machine you wish to build for +# (one of SUN, IPSC, KSR) +# + TARGET = SUN + +# +# Define SUBDIRS to be list of subdirectories of SRC to be made +# +# The include directory should be first so that the include +# files are all present and correct before any compilation +# + SUBDIRS = include global db rtdb basis inp util geom input ma tcgmsg + +# +# Define LIBPATH to be paths for libraries that you are linking in +# from precompiled sources and are not building now. These libraries +# will be searched AFTER anything you are building now. +# e.g. LIBPATH = -L/msrc/proj/mss/lib +# + LIBPATH = + +# +# Define INCPATH to be directories to get includes for +# libraries that you are not building now. These directories +# will be searched AFTER anything you are building now. +# + INCPATH = + + +########################################################## +# # +# Should NOT need to modify below here unless porting to # +# a new machine or changing compiler options # +# # +########################################################## + +# !!!! Only the SUN version is up to date !!!!! + +ifeq ($(TARGET),SUN) +# +# Sun running SunOS +# + FC = f77 + CC = gcc + AR = ar + RANLIB = ranlib + SHELL = /bin/sh + MAKE = make + MAKEFLAGS = -j 2 + INSTALL = echo $@ is built + + FOPT = -g -u -Nl99 + FOPT_REN = $(FOPT) + COPT = -g + FLDOPT = $(FOPT) + CLDOPT = $(COPT) + INCLUDES = -I. $(LIB_INCLUDES) -I$(INCDIR) $(INCPATH) + WARNINGS = -Wall +#-Wshadow -Wcast-qual -Wwrite-strings -Wpointer-arith + DEFINES = -DSUN $(LIB_DEFINES) + FFLAGS = $(FOPT) $(INCLUDES) $(DEFINES) + CFLAGS = $(COPT) $(INCLUDES) $(DEFINES) $(WARNINGS) + ARFLAGS = rcv + + LIBS = -L$(LIBDIR) $(LIBPATH) \ + -linput -lgeom -lbasis -lutil -lglobal -lrtdb -ldb -linp \ + -lutil -lma -ltcgmsg + + EXPLICITF = FALSE +endif + +ifeq ($(TARGET),IPSC) +# +# DELTA/IPSC running NX +# + FC = if77 + CC = icc + CPP = /usr/lib/cpp + AR = ar860 + + RANLIB = echo + SHELL = /bin/sh + INSTALL = rcp $@ delta2: + FOPT = -O2 -Knoieee -Mquad -node -Minline=100 + FOPT_REN = -O2 -Knoieee -Mquad -Mreentrant -Mrecursive -node + COPT = -O2 -Knoieee -Mreentrant -node + INCLUDES = -I. -I$(SRCDIR)/rtdb -I$(SRCDIR)/global -I$(SRCDIR)/tcgmsg -I$(SRCDIR)/ints \ + -I$(SRCDIR)/util -I$(SRCDIR)/ma -I$(SRCDIR)/db -I$(SRCDIR)/tcgmsg/ipcv4.0 + DEFINES = -DNX -DIPSC -DNO_BCOPY $(LIB_DEFINES) +# -DGA_TRACE + FFLAGS = $(FOPT) + CFLAGS = $(COPT) $(INCLUDES) $(DEFINES) + MAKEFLAGS = -j 2 + FLDOPT = $(FOPT) -node + CLDOPT = $(COPT) -node + ARFLAGS = rcv + LIBS = $(SRCDIR)/input/libinput.a \ + $(SRCDIR)/ddscf/libddscf.a \ + $(SRCDIR)/ints/libints.a \ + $(SRCDIR)/rtdb/librtdb.a \ + $(SRCDIR)/db/libdb.a \ + $(SRCDIR)/global/libglobal.a \ + $(SRCDIR)/trace/libtrace.a \ + $(SRCDIR)/tcgmsg/ipcv4.0/libtcgmsg.a \ + $(SRCDIR)/util/libutil.a \ + $(SRCDIR)/ma/libma.a \ + $(SRCDIR)/peigs1.0/libpeigs.a \ + $(SRCDIR)/peigs1.0/liblapack.a \ + -lkmath + + EXPLICITF = TRUE +endif + + +ifeq ($(TARGET),IBM) +# +# IBM AIX .... NOT YET TESTED !!!!! +# +# FC = xlf +# CC = xlc +# AR = ar +# RANLIB = ranlib +# INSTALL = echo +# SHELL = /bin/sh +# FOPT = -g +# COPT = -g +# INCLUDES = -I. -I../ma +# DEFINES = -DTCGMSG +# FFLAGS = -qEXTNAME $(FOPT) +# FLDOPT = $(FOPT) -b rename:.exit_,.exit +# CFLAGS = $(COPT) $(INCLUDES) $(DEFINES) +# CLDOPT = $(COPT) +# ARFLAGS = rcv +# LIBS = ../tcgmsg/ipcv4.0/libtcgmsg.a ../ma/libma.a -lc +# EXPLICITF = TRUE +# +endif + +ifeq ($(EXPLICITF),TRUE) +# +# Needed on machines where FCC does not preprocess .F files +# with CPP to get .f files +# +.SUFFIXES: +.SUFFIXES: .o .s .F .f .c + +.F.o: + $(MAKE) $*.f + $(FC) -c $(FFLAGS) $*.f + /bin/rm -f $*.f + +.F.f: + $(CPP) $(INCLUDES) $(DEFINES) < $*.F | sed '/^#/D' > $*.f + +.c.o: + $(CC) $(CFLAGS) -c $*.c +endif \ No newline at end of file diff --git a/src/config/makelib.h b/src/config/makelib.h new file mode 100644 index 0000000..098a434 --- /dev/null +++ b/src/config/makelib.h @@ -0,0 +1,59 @@ +# $Id: makelib.h,v 1.1.1.1 1994-03-29 06:44:24 d3g681 Exp $ + +# +# A makefile for a library should +# +# 1) include ../config/makefile.h ... amoung other things this will +# define TARGET from which any machine dependent actions are driven +# 2) define LIBRARY as the name of the library to be made +# 3) define OBJ as the list of object files to be made +# 4) define HEADERS as the list of header/include files to be exported +# into the common include directory +# 5) optionally define LIB_TARGETS as any additional files made in +# this subdirectory that may need cleaning up +# 6) optionally define LIB_DEFINES as any additional defines for +# the C preprocessor +# 7) optionally define LIB_INCLUDES as any additional includes +# 8) include ../config/makelib.h +# 9) define any additional targets (e.g., test programs) +# +# E.g. +# +# include ../config/makefile.h +# +# OBJ = a.o b.o c.o +# LIBRARY = libsimple.a +# HEADERS = simple.h +# LIB_TARGETS = test.o test.x +# LIB_DEFINES = -DGOODBYE="\"Have a nice day\"" +# LIB_INCLUDES = -I../testdir +# +# include ../config/makelib.h +# +# test: test.o $(LIBRARY) +# $(CC) -o $@ $^ +# +# a.o b.o c.o test.o: simple.h +# + +$(LIBRARY): $(OBJ) + /bin/rm -f $@ + $(AR) $(ARFLAGS) $@ $(OBJ) + $(RANLIB) $@ + cp -p $(LIBRARY) $(LIBDIR) + +ifdef HEADERS +include_stamp: $(HEADERS) + cp -p $(HEADERS) $(INCDIR) + touch include_stamp +else +include_stamp: + touch include_stamp +endif + +clean: + /bin/rm -f $(LIBRARY) $(OBJ) core include_stamp $(LIB_TARGETS) + + +realclean: clean + /bin/rm -f *~ \#*\# \ No newline at end of file diff --git a/src/geom/GNUmakefile b/src/geom/GNUmakefile new file mode 100644 index 0000000..589cb80 --- /dev/null +++ b/src/geom/GNUmakefile @@ -0,0 +1,9 @@ + + OBJ = geom.o geom_input.o + LIBRARY = libgeom.a + HEADERS = geom.h geomP.h + +include ../config/makefile.h +include ../config/makelib.h + +geom_input.o geom.o: geomP.h diff --git a/src/geom/geom.c b/src/geom/geom.c new file mode 100644 index 0000000..f9d5a96 --- /dev/null +++ b/src/geom/geom.c @@ -0,0 +1,683 @@ +#include +#include +#include +#include + +#include "inp.h" +#include "rtdb.h" +#include "util.h" +#include "geomP.h" +#include "tcgmsg.h" +#include "context.h" + +int ngeom_rtdb = 0; +bool active[max_geom] = {false}; +char *symbols = { + "H ", "He", "Li", "Be", "B ", "C ", "N ", "O ", "F ", "Ne", + "Na", "Mg", "Al", "Si", "P ", "S ", "Cl", "Ar", "K ", "Ca", + "Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", + "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y ", "Zr", + "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", + "Sb", "Te", "I ", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", + "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", + "Lu", "Hf", "Ta", "W ", "Re", "Os", "Ir", "Pt", "Au", "Hg", + "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", + "Pa", "U ", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", + "Md", "No", "Lr" +}; + +char *elements = { + "Hydrogen", "Helium", "Lithium", "Beryllium", "Boron", + "Carbon", "Nitrogen", "Oxygen", "Fluorine", "Neon", "Sodium", + "Magnesium", "Aluminium", "Silicon", "Phosphorous", + "Sulphur", "Chlorine", "Argon", "Potassium", "Calcium", + "Scandium", "Titanium", "Vanadium", "Chromium", "Manganese", + "Iron", "Cobalt", "Nickel", "Copper", "Zinc", "Gallium", + "Germanium", "Arsenic", "Selenium", "Bromine", "Krypton", + "Rubidium", "Strontium", "Yttrium", "Zirconium", "Niobium", + "Molybdenum", "Technetium", "Ruthenium", "Rhodium", + "Palladium", "Silver", "Cadmium", "Indium", "Tin", + "Antinomy", "Tellurium", "Iodine", "Xenon", "Caesium", + "Barium", "Lanthanum", "Cerium", "Praseodymium", "Neodymium", + "Promethium", "Samarium", "Europium", "Gadolinium", + "Terbium", "Dysprosium", "Holmium", "Erbium", "Thulium", + "Ytterbium", "Lutetium", "Hafnium", "Tantalum", "Tungsten", + "Rhenium", "Osmium", "Iridium", "Platinum", "Gold", + "Mercury", "Thallium", "Lead", "Bismuth", "Polonium", + "Astatine", "Radon", "Francium", "Radium", "Actinium", + "Thorium", "Protoactinium", "Uranium", "Neptunium", + "Plutonium", "Americium", "Curium", "Berkelium", + "Californium", "Einsteinium", "Fermium", "Mendelevium", + "Nobelium", "Lawrencium" +}; + +bool geom_check_handle(FILE *geom, char *msg) { + + bool ret_val; + + ret_val = geom > 0 && geom < max_geom; + if (ret_val) ret_val = ret_val && active[geom]; + + if (!ret_val) { + printf("%s: geometry handle invalid %d", msg, geom); + geom_err_info(msg); + } + + return ret_val; +} + +bool geom_check_cent(FILE *geom, char *msg, int icent) { + bool ret_val; + ret_val = icent > 0 && icent <= ncenter[geom]; + if (!ret_val) { + printf("%s: icent invalid %d %s\n", msg, icent, names[geom]); + geom_err_info(msg); + geom_print(geom); + } + return ret_val; +} + +bool geom_rtdb_in(FILE *rtdb) { + + /* + load in info about known geometries ... this is more + for diagnostic and debugging purposes + */ + FILE *geom; + bool ret_val = false; + int ngeom_rtdb = 0; + + if (rtdb_par_get(rtdb, "geometry:ngeom", MT_INT, 1, &ngeom_rtdb)) { + if (!rtdb_par_cget(rtdb, "geometry:names", max_geom, names_rtdb)) { + printf("geom_rtdb_in: rtdb corrupt\n"); + } else { + for (geom = 0; geom < ngeom_rtdb; geom++) { + lenr[geom] = inp_strlen(names_rtdb[geom]); + } + ret_val = true; + } + } + + return ret_val; +} + +bool geom_rtdb_out(FILE *rtdb) { + + bool ret_val; + + // output to rtdb info about known geometries + + ret_val = rtdb_par_put(rtdb, 'geometry:ngeom', MT_INT, 1, ngeom_rtdb) + && rtdb_par_cput(rtdb, 'geometry:names', max_geom, names_rtdb); + if (!ret_val) printf(" geom_rtdb_out: rtdb is corrupt "); + +} + +bool geom_rtdb_add(FILE *rtdb, char *name) { + FILE *geom; + bool status, ret_val; + int ln; + + // See if name is on the rtdb already + ln = strlen(name); + status = geom_rtdb_in(rtdb); + ret_val = true; + for (geom = 0; geom < ngeom_rtdb; geom++) { + if (strncmp(name, names_rtdb[geom], ln) == 0) { + return true; + } + } + + // Name is not present ... add and rewrite info + if (ngeom_rtdb == max_geom_rtdb) { + printf(" geom_rtdb_add: too many geometries on rtdb %s\n", name); + return false; + } + + ngeom_rtdb++; + strncpy(names_rtdb[ngeom_rtdb], name, ln); + lenr[ngeom_rtdb] = ln; + + if (!geom_rtdb_out(rtdb)) { + printf(" geom_rtdb_add: rtdb error adding %.*s\n", ln, name); + return false; + } + + return true; +} + +bool geom_err_info(char *info) { + FILE *geom; + int ngeom = 0; + /* + For internal use of the geom routines only: print out + info of known geometries to aid in diagnosing a problem + */ + for (geom = 0; geom < max_geom; geom++) { + if (active[geom]) { + ngeom++; + } + } + + printf(" %s: open geometries: %d\n", info, ngeom); + + ngeom = 0; + for (geom = 0; geom < max_geom; geom++) { + if (active[geom]) { + printf(" %d %s: \"%s\" -> \"%s\"\n", ngeom, info, names[geom], trans[geom]); + } + } + + if (ngeom_rtdb > 0) { + printf(" %s: geometries in last accessed data base: %d\n", info, ngeom_rtdb); + for (geom = 0; geom < ngeom_rtdb; geom++) { + printf(" %s\n", names_rtdb[geom]); + } + } + + return true; +} + +bool geom_rtdb_load(FILE *rtdb, FILE *geom, char *name) { + char tmp[256]; + int k; + bool status, ret_val; + + ret_val = geom_check_handle(geom, "geom_rtdb_load"); + if (!ret_val) return false; + status = geom_rtdb_in(rtdb); + + // Translate the provided name + strcpy(names[geom], name); + lenn[geom] = strlen(name); + strcpy(trans[geom], "junk"); + if (!context_rtdb_match(rtdb, name, trans[geom])) + strcpy(trans[geom], name); + lent[geom] = strlen(trans[geom]); + + // Now get the info from the data base + strcpy(tmp, "geometry:"); + strncat(tmp, trans[geom], lent[geom]); + k = strlen(tmp) + 1; + status = true; + + strcpy(tmp + k, ":ncenter"); + status = status && rtdb_par_get(rtdb, tmp, MT_INT, 1, &ncenter[geom]); + strcpy(tmp + k, ":coords"); + status = status && rtdb_par_get(rtdb, tmp, MT_DBL, max_cent * 3, &coords[1][1][geom]); + strcpy(tmp + k, ":charges"); + status = status && rtdb_par_get(rtdb, tmp, MT_DBL, max_cent, &charge[1][geom]); + strcpy(tmp + k, ":efield"); + status = status && rtdb_par_get(rtdb, tmp, MT_DBL, 3, &efield[1][geom]); + strcpy(tmp + k, ":latvec"); + status = status && rtdb_par_get(rtdb, tmp, MT_DBL, 9, &latvec[1][1][geom]); + strcpy(tmp + k, ":tags"); + status = status && rtdb_par_cget(rtdb, tmp, max_cent, &tags[1][geom]); + + if (!status) { + printf(" geom_rtdb_load: not found or rtdb corrupt: %.*s -> %.*s\n", lenn[geom], names[geom], lent[geom], trans[geom]); + geom_err_info("geom_rtdb_load"); + return false; + } + + // Determine if system is periodic or if external fields are applied + oefield[geom] = ddot(3, efield[1][geom], 1, efield[1][geom], 1) > 0.0; + operiodic[geom] = ddot(9, latvec[1][1][geom], 1, latvec[1][1][geom], 1) > 0.0; + + // Compute effective nuclear repulsion energy, dipole and interaction with external fields + geom_compute_values(geom); + + active[geom] = true; + return true; +} + +double geom_compute_values(FILE *geom) { + /* + compute effective nuclear repulsion energy, dipole and + interaction with external fields + */ + + double e, r; + int i, j; + + e = 0.0; + ndipole[0][geom] = 0.0; + ndipole[1][geom] = 0.0; + ndipole[2][geom] = 0.0; + + // compute nuclear dipole moment and usual nuclear repulsion energy + for (i = 0; i < ncenter[geom]; i++) { + for (j = 0; j < 3; j++) { + ndipole[j][geom] += charge[i][geom] * coords[j][i][geom]; + } + for (j = i + 1; j < ncenter[geom]; j++) { + r = sqrt(pow(coords[0][i][geom] - coords[0][j][geom], 2) + + pow(coords[1][i][geom] - coords[1][j][geom], 2) + + pow(coords[2][i][geom] - coords[2][j][geom], 2)); + e += charge[i][geom] * charge[j][geom] / r; + } + } + + // add in interaction of nuclear dipole with external field + e += ddot(3, ndipole[0][geom], 1, efield[0][geom], 1); + + erep[geom] = e; +} + +bool geom_rtdb_store(FILE *rtdb, char *name, FILE *geom) { + bool status, ret_val; + char tmp[256]; + + ret_val = geom_check_handle(geom, "geom_rtdb_store"); + if (!ret_val) { + return false; + } + + // Update the name if provided + if (name != NULL && strcmp(name, "") != 0) { + strcpy(names[geom], name); + lenn[geom] = strlen(name); + } + + // If not process 0 return ... this is so that input routines can be completely single threaded + if (nodeid() != 0) { + return true; + } + + // Try to translate the name + strcpy(trans[geom], "junk"); + if (!context_rtdb_match(rtdb, name, trans[geom])) { + strcpy(trans[geom], name); + } + lent[geom] = strlen(trans[geom]); + + // Now put the info into the data base + strcpy(tmp, "geometry:"); + k = strlen(tmp); + status = true; + + strcpy(tmp[k], ":ncenter\0"); + status = status && rtdb_par_put(rtdb, tmp, MT_INT, 1, &ncenter[geom]); + strcpy(tmp[k], ":coords\0"); + status = status && rtdb_par_put(rtdb, tmp, MT_DBL, ncenter[geom] * 3, &coords[1][1][geom]); + strcpy(tmp[k], ":charges\0"); + status = status && rtdb_par_put(rtdb, tmp, MT_DBL, ncenter[geom], &charge[1][geom]); + strcpy(tmp[k], ":efield\0"); + status = status && rtdb_par_put(rtdb, tmp, MT_DBL, 3, &efield[1][geom]); + strcpy(tmp[k], ":latvec\0"); + status = status && rtdb_par_put(rtdb, tmp, MT_DBL, 9, &latvec[1][1][geom]); + strcpy(tmp[k], ":tags\0"); + status = status && rtdb_par_cput(rtdb, tmp, ncenter[geom], &tags[1][geom]); + + // Insert translated name into list of known geometries + status = status && geom_rtdb_add(rtdb, name); + + // Check that all rtdb operations were successful + if (!status) { + printf(" geom_rtdb_store: write to rtdb failed: %.*s -> %.*s\n", lenn[geom], names[geom], lent[geom], trans[geom]); + geom_err_info("geom_rtdb_store"); + return false; + } + + return true; +} + +bool geom_rtdb_delete(FILE *rtdb, char *name) { + char translation[256], tmp[256]; + int lt, geom, geom2, k; + bool status, set_true = false; + + // try to translate the provided name + if (rtdb_par_cget(rtdb, name, 1, translation) != 0) { + strcpy(translation, name); + } + lt = strlen(translation); + + // locate name in list and remove + status = geom_rtdb_in(rtdb); + for (geom = 1; geom <= ngeom_rtdb; geom++) { + if (strncmp(names_rtdb[geom], translation, lt) == 0) { + set_true = true; + break; + } + } + + if (!set_true) { + printf(" geom_rtdb_delete: no such geometry %.*s -> %.*s\n", strlen(name), name, lt, translation); + } + + for (geom2 = geom + 1; geom2 <= ngeom_rtdb; geom2++) { + strcpy(names_rtdb[geom2 - 1], names_rtdb[geom2]); + } + ngeom_rtdb--; + + status = geom_rtdb_out(rtdb); + + // delete each entry associated with a geometry in the database + strcpy(tmp, "geometry:"); + k = strlen(tmp); + + strcpy(tmp[k], ":ncenter"); + status = status && rtdb_par_delete(rtdb, tmp); + strcpy(tmp[k], ":coords"); + status = status && rtdb_par_delete(rtdb, tmp); + strcpy(tmp[k], ":charges"); + status = status && rtdb_par_delete(rtdb, tmp); + strcpy(tmp[k], ":efield"); + status = status && rtdb_par_delete(rtdb, tmp); + strcpy(tmp[k], ":latvec"); + status = status && rtdb_par_delete(rtdb, tmp); + strcpy(tmp[k], ":tags"); + status = status && rtdb_par_delete(rtdb, tmp); + + // check status of all rtdb stores + if (!status) { + printf(" geom_rtdb_delete: rtdb corrupt %.*s\n", lt, translation); + geom_err_info("geom_rtdb_delete"); + return false; + } + + return true; +} + +bool geom_create(FILE *geom, char *name) { + char translation[256]; + int i; + + // Assign the next free slot for a geometry + for (i = 0; i < max_geom; i++) { + if (!active[i]) { + break; + } + } + + if (i == max_geom) { + printf("geom_create: too many geoms trying to create %s\n", name); + geom_err_info("geom_create"); + return false; + } + + // Store info about the geometry + strcpy(names[i], name); + strcpy(trans[i], " "); + lenn[i] = strlen(name); + ncenter[i] = 0; + active[i] = true; + + return true; +} + +bool geom_destroy(FILE *geom) { + + bool ret_val; + bool ret_val = geom_check_handle(geom, "geom_destroy"); + if (!ret_val) return false; + + active[geom] = false; + + return true; +} + +bool geom_cart_set(FILE *geom, int ncent, char *t[], double c[3][ncent], double q[ncent]) { + + int i; + bool ret_val; + + if (!geom_check_handle(geom, "geom_cart_set")) { + return false; + } + + if (ncent <= 0 || ncent > max_cent) { + printf("geom_cart_set: too many centers %d %s\n", ncent, names[geom]); + return false; + } + + ncenter[geom] = ncent; + for (i = 0; i < ncent; i++) { + strcpy(tags[i][geom], t[i]); + charge[i][geom] = q[i]; + coords[0][i][geom] = c[0][i]; + coords[1][i][geom] = c[1][i]; + coords[2][i][geom] = c[2][i]; + } + /* + compute effective nuclear repulsion energy, dipole and + interaction with external fields + */ + geom_compute_values(geom); + + return true; +} + +bool geom_cart_get(FILE *geom, int *ncent, char *t[], double c[][ncent], double q[]) { + + int i; + + if (!geom_check_handle(geom, "geom_cart_get")) { + return false; + } + + *ncent = ncenter[geom]; + for (i = 0; i < *ncent; i++) { + strcpy(t[i], tags[i][geom]); + q[i] = charge[i][geom]; + c[0][i] = coords[0][i][geom]; + c[1][i] = coords[1][i][geom]; + c[2][i] = coords[2][i][geom]; + } + + return true; +} + +bool geom_cent_get(FILE *geom, int icent, char *t, double c[3], double *q) { + bool ret_val; + ret_val = geom_check_handle(geom, "geom_cent_get"); + if (!ret_val) return false; + ret_val = geom_check_cent(geom, "geom_cent_get", icent); + if (!ret_val) return false; + + strcpy(t, tags[icent][geom]); + c[0] = coords[0][icent][geom]; + c[1] = coords[1][icent][geom]; + c[2] = coords[2][icent][geom]; + q = charge[icent][geom]; + + return true; +} + +bool geom_cent_set(FILE *geom, int icent, char *t, double c[3], double q) { + if (!geom_check_handle(geom, "geom_cent_set")) { + return false; + } + + if (!geom_check_cent(geom, "geom_cent_set", icent)) { + return false; + } + + strcpy(tags[icent][geom], t); + coords[0][icent][geom] = c[0]; + coords[1][icent][geom] = c[1]; + coords[2][icent][geom] = c[2]; + charge[icent][geom] = q; + + geom_compute_values(geom); + + return true; +} + +bool geom_ncent(FILE *geom, int ncent) { + bool ret_val; + ret_val = geom_check_handle(geom, "geom_ncent"); + if (!ret_val) return false; + ncent = ncenter[geom]; + + return true; +} + +bool geom_cent_tag(FILE *geom, int icent, char *tag) { + if (!geom_check_handle(geom, "geom_cent_tag")) { + return false; + } + + if (!geom_check_cent(geom, "geom_cent_tag", icent)) { + return false; + } + + strcpy(tag, tags[icent][geom]); + + return true; +} + +bool geom_latvec_set(FILE *geom, double vectors[3][3]) { + errquit("geom_latvec_set: not yet!", 0); + + return false; +} + +bool geom_latvec_get(FILE *geom, double vectors[3][3]) { + errquit("geom_latvec_get: not yet!", 0); + + return false; +} + +bool geom_efield_set(FILE *geom, double efield[3]) { + errquit("geom_efield_set: not yet!", 0); + geom_set_values(geom); + + return false; +} + +bool geom_efield_get(FILE *geom, double efield) { + errquit("geom_efield_get: not yet!", 0); + + return false; + +} + +bool geom_print(FILE *geom) { + /* + Basic printing of cartesian geometry ... needs support for + user defined units, internal coords, different formats, ... + */ + + int icent, ivec, i; + bool ret_val; + + if (!geom_check_handle(geom, "geom_print")) { + return false; + } + + printf(" Geometry (au) \"%s\" -> \"%s\"\n", names[geom], trans[geom]); + printf(" -------------\n"); + printf(" No. Tag Charge X Y Z\n"); + printf(" ---- ---------------- ---------- -------------- -------------- --------------\n"); + + for (icent = 0; icent < ncenter[geom]; icent++) { + printf("%4d %16s %10.6f %14.8f %14.8f %14.8f\n", icent, tags[icent][geom], charge[icent][geom], coords[0][icent][geom], coords[1][icent][geom], coords[2][icent][geom]); + } + + printf("Effective nucler repulsion charge (au) %18.10f", erep[geom]); + + if (operiodic[geom]) { + printf("Periodic lattice vectors (au)\n"); + printf(" -----------------------------\n"); + printf(" X Y Z\n"); + printf(" ---------------- ---------------- ----------------\n"); + for (ivec = 0; ivec < 3; ivec++) { + printf(" %17.10f %17.10f %17.10f\n", latvec[ivec][0][geom], latvec[ivec][1][geom], latvec[ivec][2][geom]); + } + + } + + if (oefield[geom]) { + printf("Electric Field (au)\n"); + printf(" -------------------\n"); + printf(" X Y Z\n"); + printf(" ---------------- ---------------- ----------------\n"); + printf(" %17.10f %17.10f %17.10f\n", efield[0][geom], efield[1][geom], efield[2][geom]); + } + + return true; +} + +bool geom_tag_to_element(char *tag, char *symbol, char *element, int *atn) { + /* + attempt to figure out which element a tag refers to + and return the symbol, name and atomic no. + */ + + bool ret_val; + int lbuf, ind; + char buf[17]; + char sym1[15] = {'h', 'b', 'c', 'n', 'o', 'f', 'p', 's', 'k', 'v', 'y', 'i', 'w', 'u'}; + int atn1[14] = {1, 5, 6, 7, 8, 9, 15, 16, 19, 23, 39, 53, 74, 92}; + + ret_val = false; + /* + eliminate conventions that refer to centers used for + computation purposes .. just bq for now + */ + lbuf = strlen(buf); + if (lbuf == 0) return false; + + for (int i = 0; i < lbuf; i++) { + buf[i] = tolower(buf[i]); + } + + if (strncmp(buf, "bq", 2) == 0) { + strcpy(element, "point charge"); + strcpy(symbol, "bq"); + atn = 0; + return false; + } + /* + Attempt to match the first 4 characters of the + full names of the elements + */ + atn = 0; + if (lbuf >= 4) { + for (int i = 0; i < nelements; i++) { + if (strncmp(buf, elements[i], 4) == 0) { + strcpy(symbol, symbols[i]); + strcpy(element, elements[i]); + atn = i; + ret_val = true; + return true; + } + } + } + /* + Failed ... attempt to match the first two characters + against two character element names + */ + if (buf[1] != ' ') { + for (int i = 0; i < nelements; i++) { + if (strncmp(buf, symbols[i], 2) == 0) { + strcpy(symbol, symbols[i]); + strcpy(element, elements[i]); + atn = i; + ret_val = true; + return true; + } + } + } + + // Last ditch attempt ... match against 1 character symbols + for (int i = 0; i < 14; i++) { + if (buf[0] == sym1[i]) { + ind = atn1[i]; + strcpy(symbol, symbols[ind]); + strcpy(element, elements[ind]); + atn = ind; + ret_val = true; + return true; + } + } + + // Nothing matched + strcpy(symbol, " "); + strcpy(element, " "); + atn = 0; + return false; + +} \ No newline at end of file diff --git a/src/geom/geom.doc b/src/geom/geom.doc new file mode 100644 index 0000000..ab12b0b --- /dev/null +++ b/src/geom/geom.doc @@ -0,0 +1,113 @@ + +The geometry data includes + + 1) A description of the coordinates of all types of centers (e.g., + atom, charge, basis function) + + 2) Charges (and I guess possibly other potentials) associated with + those centers + + 3) Tags (names) of centers + + 4) Masses associated with centers + + 5) Variables for optimization (e.g., via constrained cartesians + or zmatrix variables) + + 6) Any other simple scalar/vector attributed associated + specifically with a center + +Operations + + 1) Store/retrieve from the database + + logical geom_rtdb_load(rtdb, name, geom) + integer rtdb [input] + character*(*) name [input] + integer geom [output] + + logical geom_rtdb_store(rtdb, 'geometry', geom) + integer rtdb [input] + character*(*) name [input] + integer geom [input] + + 2) Create/destroy + + logical geom_create(geom) + integer geom [output] + + logical geom_destroy(geom) + integer geom [input] + + 3) Set/get commmon values for all centers + + logical geom_cart_set(geom, ncent, tags, coords, charges) + integer geom [input] + integer ncent [input] + character*(*) tags(ncent) [input] + character*(*) coords(3, ncent) [input] + character*(*) charges(ncent) [input] + + logical geom_cart_get(geom, ncent, tags, coords, charges) + integer geom [input] + integer ncent [output] + character*(*) tags(ncent) [output] + character*(*) coords(3, ncent) [output] + character*(*) charges(ncent) [output] + + 4) Set/get common values for specific centers + + logical geom_cent_set(geom, icent, tag, coord, charge) + integer geom [input] + integer ncent [input] + character*(*) tag [input] + character*(*) coords(3) [input] + character*(*) charge [input] + + logical geom_cent_get(geom, icent, tag, coord, charge) + integer geom [input] + integer ncent [output] + character*(*) tag [output] + character*(*) coords(3) [output] + character*(*) charge [output] + + 5) Inquiry routines + + integer function geom_ncent(geom) + integer geom [input] + + logical function geom_cent_tag(geom, icent, tag) + integer geom [input] + integer icent [input] + character*(*) tag [output] + + 6) Set/get specific values for specific centers + + There are two possibilities here + + a) adopt an extensible definition of properties associated + with a center. This includes registering new properties + with a name and routines to set/get/load/store the values + and some general format (e.g., netcdf) for describing + and passing data. + + b) adopt a static definition of the data structures and + require recompilation after the structures have been changed + and new routines provided. + + Do we think that new properties will be added very regularly or + that this will become very infrequent? I tend to think the latter, + so a) is not yet worth the effort. Since b) requires very little + effort we can always change our minds and do a) later. + + + 7) Zmatrix routines ... not yet defined + + n_zmat_cent, n_zmat_vars, ... + logical geom_zmat_defined() + call geom_zmat_get + call geom_zmat_set + ... + + +Data on the rtdb diff --git a/src/geom/geom.h b/src/geom/geom.h new file mode 100644 index 0000000..1da34ee --- /dev/null +++ b/src/geom/geom.h @@ -0,0 +1,31 @@ +#ifndef _GEOM_H +#define _GEOM_H + +#include + + bool geom_check_handle(FILE *, char *); + bool geom_check_cent(FILE *, char *, int); + bool geom_rtdb_in(FILE *); + bool geom_rtdb_out(FILE *); + bool geom_rtdb_add(FILE *, char *); + bool geom_err_info(char *); + bool geom_rtdb_load(FILE *, FILE *, char *); + double geom_compute_values(FILE *); + bool geom_rtdb_store(FILE *, char *, FILE *); + bool geom_rtdb_delete(FILE *, char *); + bool geom_create(FILE *, char *); + bool geom_destroy(FILE *); + bool geom_cart_set(FILE *, int, char *[], double[3][], double[]); + bool geom_cart_get(FILE *, int *, char *[], double *[], double []); + bool geom_cent_get(FILE *, int, char *, double[3], double *); + bool geom_cent_set(FILE *, int, char *, double[3], double); + bool geom_ncent(FILE *, int); + bool geom_cent_tag(FILE *, int, char *); + bool geom_latvec_set(FILE *, double [3][3]); + bool geom_latvec_get(FILE *, double [3][3]); + bool geom_efield_set(FILE *, double[3]); + bool geom_efield_get(FILE *, double); + bool geom_print(FILE *); + bool geom_tag_to_element(char *, char *, char *, int *); + +#endif \ No newline at end of file diff --git a/src/geom/geomP.h b/src/geom/geomP.h new file mode 100644 index 0000000..c4aee15 --- /dev/null +++ b/src/geom/geomP.h @@ -0,0 +1,84 @@ +#ifndef _GEOMP_H +#define _GEOMP_H +/* + Private fortran include file for the geometry routines + + Parameters + + max_geom = maximum no. of geometries + max_cent = maximum no. of centers in a geometry + max_geom_rtdb = maximum no. of geometries stored in the rtdb + nelments = no. of elements that info is stored about + + [The only thing that cannot be dynamically allocated are the + character variables for the tags ... I was lazy and just statically + dimensioned everything ... just drudge work to dynamically + allocate though if necessary ... which it hopefully won't be + ... since only geom.F (and maybe the basis routines) include + this header file only these need to be recompiled if the parameters + are changed] + + Members of /cgeometry/ + + ngeom_rtdb = current no. of geometries on the rtdb + active(1:max_geom) = true if this geometry is open + ncenter(1:max_geom) = no. of centers in this geometry + coords(1:3,1:max_cent,1:max_geom) = cartesian coords of this geometry + charge(1:max_cent,1:max_geom) = charges associated with centers + dipole ... not yet + quadrupole ... not yet + pseudopotential ... not yet + efield(1:3,1:max_geom) = external electric field applied to this system + oefield = true if efield is on + latvec(1:3,1:3,1:max_geom) = vectors specifing periodicity + (null vector gives no periodicity) + operiodic = true if a lattice vector is non-null + erep(1:max_geom) = interaction energy of centers with each other + and external fields. At its simplest this is + just the nuclear repulsion energy + ndipole(1:3,1:max_geom) = nuclear dipole moment + Members of /ccgeometry/ + + names(1:max_geom) = names of open geometries + trans(1:max_geom) = translations of names of open geoms + names_rtdb(1:max_geom) = names of geometries in the rtdb + tag(1:max_cent,1:max_geom) = tags associated with centers + lenn(1:max_geom) = length of names(geom) minus trailing blanks + lent(1:max_geom) = length of trans(geom) ... + lenr(1:max_geom) = length of names_rtdb(geom) ... + symbols(1:nelements) = symbols for elements + elements(1:nelements) = names of elements +*/ + +#define MAX_GEOM 2 +#define MAX_CENT 1000 +#define MAX_GEOM_RTDB 100 +#define NELEMENTS 103 + +typedef struct { + double coords[3][MAX_CENT][MAX_GEOM]; + double charge[MAX_CENT][MAX_GEOM]; + double efield[3][MAX_GEOM]; + double latvec[3][3][MAX_GEOM]; + double erep[MAX_GEOM]; + double ndipole[3][MAX_GEOM]; + int ncenter[MAX_GEOM]; + int active[MAX_GEOM]; + int lenn[MAX_GEOM]; + int lent[MAX_GEOM]; + int lenr[MAX_GEOM]; + int operiodic[MAX_GEOM]; + int oefield[MAX_GEOM]; + int ngeom_rtdb; +} CGeometry; + +typedef struct { + char names[MAX_GEOM][256]; + char trans[MAX_GEOM][256]; + char names_rtdb[MAX_GEOM_RTDB][256]; + char tags[MAX_CENT][MAX_GEOM][16]; + char symbols[NELEMENTS][2]; + char elements[NELEMENTS][16]; +} CCGeometry; + +#endif // _GEOMP_H diff --git a/src/geom/geom_input.c b/src/geom/geom_input.c new file mode 100644 index 0000000..93ef968 --- /dev/null +++ b/src/geom/geom_input.c @@ -0,0 +1,98 @@ +#include +#include + +#include "inp.h" +#include "geom.h" +#include "tcgmsg.h" + +void geom_input(FILE *rtdb, bool print) { + char field[255]; // for character input + char name[255]; // for name of geometry + char units[12]; // holds units of coordinates + int ncenter; // counts no. of centers as input + FILE *geom; // handle for geometry + bool status; // scratch for return codes + const int max_center = 1000; // parameter for local array dimension + double coords[max_center][3]; + double charge[max_center]; + char tags[max_center][16]; + /* + read a geometry from the input deck and output it + to the rtdb. + + current input line should begin 'geometry ...' + + Cartesians only for now + */ + if (nodeid() != 0) return; + + // Check that this is indeed a geometry line + inp_set_field(0); // goto start of line + if (!inp_a(field)) { + errquit("geom_input: no input present", 0); + } + if (!inp_compare(false, "geom", field)) { + errquit("geom_input: not geometry input", 0); + } + + // geometry [] [units ] + strcpy(units, "atomic units"); + strcpy(name, " "); + while (inp_a(field)) { + if (inp_compare(false, "units", field)) { + if (!inp_a(units)) { + errquit("geom_input: geometry [] [units ]", 0); + } + geom_check_units(units); + } else { + if (strcmp(name, " ") != 0) { + errquit("geom_input: geometry [] [units ]", 0); + } + strcpy(name, field); + } + } + + if (!geom_create(geom, name)) { + errquit("geom_input: geom_create failed !", 0); + } + + // tag charge x y z + ncenter = 0; + while (inp_read()) { + status = inp_a(field); + if (inp_compare(false, "end", field)) break; + + else { + if ((ncenter + 1) == max_center) + errquit("geom_input: too many centers?", ncenter); + strcpy(tags[ncenter + 1], field); + status = status && inp_f(charge[ncenter + 1]); + status = status && inp_f(coords[0][ncenter + 1]); + status = status && inp_f(coords[1][ncenter + 1]); + status = status && inp_f(coords[2][ncenter + 1]); + if (!status) + errquit("geom_input: ", 0); + ncenter++; + } + } + + if (!geom_cart_set(geom, ncenter, tags, coords, charge)) { + errquit("geom_input: geom_cart_set failed", 0); + } + + if (print) { + if (!geom_print(geom)){ + errquit("geom_input: print failed ", 0); + } + } + + if (!geom_rtdb_store(rtdb, name, geom)) { + errquit("geom_input: geom_rtdb_store failed", 0); + } + + if (!geom_destroy(geom)) { + errquit("geom_input: geom_destroy failed", 0); + } + + // done +} diff --git a/src/include/GNUmakefile b/src/include/GNUmakefile new file mode 100644 index 0000000..0548e7d --- /dev/null +++ b/src/include/GNUmakefile @@ -0,0 +1,20 @@ +# +# This directory is a central repository for all include +# files. The makefile in each subdirectory should contain +# a rule that keeps this directory up to date +# + +include ../config/makefile.h + +includes: + for dir in $(SUBDIRS); do \ + echo Making include_stamp in $(SRCDIR)/$$dir ; \ + (cd $(SRCDIR)/$$dir; $(MAKE) include_stamp) ; \ + done + +include_stamp: + echo Nothing to be done + +realclean clean: + echo Header files not removed + /bin/rm -f *~ \#*\# \ No newline at end of file diff --git a/src/inp/GNUmakefile b/src/inp/GNUmakefile new file mode 100644 index 0000000..94c090a --- /dev/null +++ b/src/inp/GNUmakefile @@ -0,0 +1,13 @@ +# $Id: GNUmakefile,v 1.1.1.1 1994-03-29 06:44:37 d3g681 Exp $ + + OBJ = inp.o + LIBRARY = libinp.a + LIB_TARGETS = test.o test + HEADERS = inp.fh + + +include ../config/makefile.h +include ../config/makelib.h + +test: test.o input.o + $(FC) $(FFLAGS) -o $@ $^ $(LIBS) diff --git a/src/inp/inp.doc b/src/inp/inp.doc new file mode 100644 index 0000000..0062062 --- /dev/null +++ b/src/inp/inp.doc @@ -0,0 +1,191 @@ + All routines are declared in the header file 'inp.h' + + + subroutine inp_init(ir, iw) + + Initialize free format input routines to take input from + fortran unit ir and send their output to fortran unit iw. + The input file is processed from the current location. + + inp_init() shuld be invoked each time the input file is + repositioned using other than inp_*() routines (e.g., rewind). + + + logical function inp_read() + + Read a line from the input and split it into white space (blank + or tab) separated fields. White space may be incorporated into a + field by enclosing it in quotes ("). The case of input is + preserved. Blank lines are ignored, and text from a pound or + hash symbol (#) to the end of the line is treated as a comment. + A backslash(\) at the end of a line (only white space may appear + after it) may be used to concatentate physical input lines into + one logical input line. A semicolon (;) may be used to split a + physical input line into multiple logical input lines. The + special meaning of hash (#), semicolon (;) and quotation (") + characters may be avoided only by prefacing them with a backslash + (this must be done even if the character is inside a quoted + character string). + + The no. of fields read is set to 0, there being a total of + inp_nfield() fields in the line. + + If a non-blank line is successfully parsed then .true. is returned. + + Otherwise an internal error message is set and .false. is returned. + + Possible errors include detection of EOF (inp_eof() may be used + to check for this condition) or failure to parse the line (e.g., + a character string without a terminating quote). + + EOF may be indicated by end of the physical input file, or by a + physical input line that begins with either asterisk (*), period + (.) or EOF (ignoring case), and has only trailing white space. + + There is a maximum input line width of 256 characters. + + + logical function inp_i(integer i) + + Attempt to read the next field as an integer. + Upon success return .true. and advance to the next field. + Otherwise return .false., save an internal error message and do + not change the current field. + + + logical function inp_f(double precision d) + + Attempt to read the next field as a floating point number. + Upon success return .true. and advance to the next field. + Otherwise return .false., save an internal error message and do + not change the current field. + + + logical function inp_a(character*(*) a) + + Attempt to read the next field as a character string. + Upon success return .true. and advance to the next field. + Otherwise return .false., save an internal error message and do + not change the current field. + + + logical function inp_a_trunc(character*(*) a) + + Attempt to read the next field as a character string, quietly + discarding any data that does not fit in the user provided buffer. + Upon success return .true. and advance to the next field. + Otherwise return .false., save an internal error message and do + not change the current field. + + logical function inp_line(character*(*) z) + character*(*) z + + Return in z as much of the entire input line as it will hold and + quietly discard any overflow. Upon success return .true., + otherwise save an internal error message and return .false. + + integer function inp_n_field() + + Returns the no. of fields in the current input line (1, ...). A + value of 0 implies either that EOF or some other error was + detected or inp_read() has not yet been called. + + integer function inp_cur_field() + + Returns the no. of fields in the input line that have been + processed so far (0, ...). Thus if inp_cur_field() returns 2, + then the next field read by inp_f() etc. will be field 3. + + subroutine inp_set_field(value) + integer value + + Sets the current field (as returned by inp_cur_field) to be + value. 0 <= value <= inp_n_field(). An out of range value + results in error termination. + + + subroutine inp_prev_field() + + A convenience routine that positions you to read the field (on + the current input line) that was last read. It is simply + implemented as + + call inp_set_field(max(0,inp_cur_field()-1)) + + At the beginning of the line this is a null operation. + + + logical function inp_compare(ocase, a, b) + logical ocase + character*(*) a, b + + Return .true. iff all the characters in A match the first + len(A) characters of B. If ocase is .true. then comparisons are + case sensitive, otherwise comparisons ignore case. + + + logical function inp_match(nrec, ocase, test, array, ind) + integer nrec + logical ocase + character*(*) test + character*(*) array(nrec) + integer ind + + Let L be the length of the character string test ignoring + trailing blanks. Attempt to find a unique match of test(1:L) + against elements of array(*). If ocase is .true. then + comparisons are case sensitive, otherwise comparisons ignore + case. + + If a unique match is made return the index of the element in ind + and return .true. + + If the match is ambiguous set ind to 0, and return .false. + + If no match is found set ind to -1 and return .false. + + + logical function inp_search(ocase, z) + character*(*) z + logical ocase + + Position the input file at the next logical input line which has + a first input field that matches the leading non-blank characters + in z. If ocase is .true. then matches are case sensitive. + + If such a line is found then return .true., and reset the current + input field to 0 (i.e., as if inp_read() had just been called). + + If no such line is found return .false.. The file will be either + at EOF or at a line which was not successfully parsed. EOF may + be detected by inp_eof(). + + + logical function inp_eof() + + Return .true. if EOF has been detected, .false. otherwise. + + + subroutine inp_lcase(z) + character*(*) z + + Lowercase the character string z + + + integer function inp_strlen(z) + character*(*) z + + Return the index of the last non-blank character in z, 0 being + returned for a fully blank string. + + + subroutine inp_errout() + + If there is an internal error message, print out its value, + the current line number and its contents. If appropriate + indicate the problematic position in the current input line. + + + subroutine inp_outrec() + + Print out the current input line. \ No newline at end of file diff --git a/src/inp/inp.h b/src/inp/inp.h new file mode 100644 index 0000000..d86cb78 --- /dev/null +++ b/src/inp/inp.h @@ -0,0 +1,18 @@ +#ifndef _INP_H +#define _INP_H + + int inp_i(); + float inp_f(); + char *inp_a(); + int inp_n_field(); + int inp_cur_field(); + int inp_match(); + int inp_read(); + int inp_search(); + int inp_compare(); + int inp_strlen(); + int inp_line(); + int inp_eof(); + int inp_a_trunc(); + +#endif \ No newline at end of file diff --git a/src/inp/inpP.h b/src/inp/inpP.h new file mode 100644 index 0000000..0ddece8 --- /dev/null +++ b/src/inp/inpP.h @@ -0,0 +1,34 @@ +#ifndef _INPP_H +#define _INPP_H + +// Private header file for free format input routines + +#include + + #define MAX_WIDTH 256 // Maximum no. of characters in an input line + #define MAX_FIELD MAX_WIDTH/2 + 1 // Maximum no. of fields in an input line + char ja[256]; // Input buffers ... MUST match max_width + char tmp[256]; // Same size work space + char errmsg[80]; // Error message + char xcomm; // Comment character + char xsplit; // Character to split physical input lines + char xback; // Backslash for concatenation and quoting + char xquote; // Quotation marks for strings + char xblnk; // Space + char xtab; // Tab + + int jrec; // No. of current field + int jump; // No. of fields in current line + int istrt[MAX_FIELD]; // Start of fields + int inumb[MAX_FIELD]; // Length of fields + int nstart[MAX_FIELD]; // Start of fields + int nend[MAX_FIELD]; // End of fields + int iwidth; // Length of current logical input line + int nline; // Current logical line inside physical line + int noline; // No. of logical lines inside physical line + int input_line; // No. of current physical input line + int nerr; // ???? + bool oswit; // True if EOF has beeen detected + int ierrpos; // Input char position where error was detected + +#endif // _INPP_H diff --git a/src/inp/test.c b/src/inp/test.c new file mode 100644 index 0000000..cc2020c --- /dev/null +++ b/src/inp/test.c @@ -0,0 +1,36 @@ +#include + +#include "inp.h" + +void test() { + char aval[30]; + int i, ival, line; + double dval; + + line = 0; + inp_init(5, 6); + + while (inp_read()) { + line = line + 1; + inp_outrec(); + for (i = 1; i <= inp_n_field(); i++) { + if (inp_i(&ival)) { + printf("line=%d, field=%d, integer=%d\n", line, i, ival); + } else if (inp_f(&dval)) { + printf("line=%d, field=%d, double=%.2lf\n", line, i, dval); + } else if (inp_a(aval)) { + printf("line=%d, field=%d, string=%s\n", line, i, aval); + } else { + printf("line=%d, field=%d, error!\n", line, i); + inp_errout(); + } + } + } + + if (inp_eof()) { + printf("EOF detected at line %d\n", line); + } else { + printf("input failed at line %d\n", line); + } + inp_errout(); +} diff --git a/src/rtdb/GNUmakefile b/src/rtdb/GNUmakefile new file mode 100644 index 0000000..084eacd --- /dev/null +++ b/src/rtdb/GNUmakefile @@ -0,0 +1,47 @@ +# $Id: GNUmakefile,v 1.1.1.1 1994-03-29 06:44:27 d3g681 Exp $ + +include ../config/makefile.h + + OBJ = rtdb_f2c.o rtdb.o rtdb_par.o rtdb_par_f2c.o \ + context.o context_f2c.o + LIBRARY = librtdb.a + + HEADERS = context.h rtdb.h rtdb.h context.h + + LIB_TARGETS = test test.o rtdbtest rtdbtest.o interact context davetest \ + context.o davetest.o interact.o rtdb_par_f2c.o \ + cntx.o cntx + + LIB_INCLUDES = -I../db + + TEST_LIBS = $(LIBRARY) $(LIBS) + +include ../config/makelib.h + +davetest: davetest.o $(LIBRARY) + $(FC) $(FFLAGS) -o $@ davetest.o librtdb.a ../db/libdb.a ../util/libutil.a ../ma/libma.a + +rtdb_par_f2c.c: rtdb_f2c.c + sed 's/rtdb_/rtdb_par_/g' rtdb_f2c.c > rtdb_par_f2c.c + +cntx: cntx.o $(LIBRARY) + $(FC) $(FFLAGS) -o $@ cntx.o $(LIBS) + +interact: interact.o $(LIBRARY) + $(CC) $(CFLAGS) -o $@ interact.o $(TEST_LIBS) + +test: test.o $(LIBRARY) + $(FC) $(FFLAGS) -o $@ test.o $(TEST_LIBS) + +rtdbtest: rtdbtest.o $(LIBRARY) + $(CC) $(CFLAGS) -o $@ rtdbtest.o $(TEST_LIBS) + +rtdbpartest: rtdb_par_test.o $(LIBRARY) + $(CC) $(CFLAGS) -o $@ rtdb_par_test.o $(TEST_LIBS) + +rtdb.o rtdbf2c.o rtdbtest.o: misc.h +rtdb.o rtdbf2c.o rtdbtest.o: rtdb.h +test.o: rtdb.h + +context: context.o $(LIBRARY) + $(CC) $(CFLAGS) -o $@ context.o $(TEST_LIBS) diff --git a/src/rtdb/context.h b/src/rtdb/context.h new file mode 100644 index 0000000..0e5deda --- /dev/null +++ b/src/rtdb/context.h @@ -0,0 +1,15 @@ +#ifndef _CONTEXT_H +#define _CONTEXT_H + +#include + +bool context_set(const char *); +char *context_get(void); +bool context_rtdb_store(int); +bool context_rtdb_load(int); +bool context_push(const char *); +bool context_pop(const char *); +bool context_rtdb_match(int, const char *, int, char *); +bool context_prefix(const char *, char *, int); + +#endif \ No newline at end of file diff --git a/src/rtdb/rtdb.h b/src/rtdb/rtdb.h new file mode 100644 index 0000000..fd015aa --- /dev/null +++ b/src/rtdb/rtdb.h @@ -0,0 +1,154 @@ +#ifndef _RTDB_H +#define _RTDB_H + +/* + All routines return TRUE (1) on success, FALSE (0) on failure. + + int rtdb_open(const char *filename, const char *mode, int *handle) + + Filename = path to file associated with the data base + mode = 'new' Open only if it does not exist already + 'old', Open only if it does exist already + 'unknown' Create new or open existing (preserving contents) + 'empty' Create new or open existing (deleting contents) + 'scratch' Create new or open existing (deleting contents) + and automatically delete upon closing. Also, items + cached in memory are not written to disk. + + handle = returns handle by which all future references to the + data base are made + + + + int rtdb_close(const int handle, const char *mode) + + Close the data base + + handle = handle to RTDB + mode = 'keep' Preserve the data base file to enable restart + 'delete' Delete the data base file freeing all resources + + mode is overridden by opening the data base with + mode='scratch' in which instance it is always deleted + upon closing + + + int rtdb_get_info(const int handle, const char *name, int *ma_type, + int *nelem, char date[26]) + + Get info about an entry from the data base + + handle = handle to RTDB + name = entry name (null terminated character string) + ma_type = returns MA type of the entry + nelem = returns no. of elements of the given type + date = returns date of insertion (null terminated character string) + + + int rtdb_put(const int handle, const char *name, const int ma_type, + const int nelem, const void *array) + + Insert an entry into the data base replacing previous entry + + handle = handle to RTDB + name = entry name (null terminated character string) + ma_type = MA type of the entry + nelem = no. of elements of the given type + array = data to be inserted + + + int rtdb_get(const int handle, const char *name, const int ma_type, + const int nelem, void *array) + + Get an entry from the data base + + handle = handle to RTDB + name = entry name (null terminated character string) + ma_type = MA type of the entry which must match entry type + nelem = size of array in units of ma_type + array = user provided buffer that returns data + + int rtdb_ma_get(const int handle, const char *name, int *ma_type, + int *nelem, int *ma_handle) + + Get an entry from the data base returning an MA handle + + handle = handle to RTDB + name = entry name (null terminated character string) + ma_type = returns MA type of the entry + nelem = returns no. of elements of type ma_type in data + ma_handle= returns MA handle to data + + int rtdb_first(const int handle, const int namelen, char *name) + + Return the name of the first (user inserted) entry in the data base. + The order is effectively random. + + handle = handle to RTDB + namelen = size of user provided buffer name + name = name of entry is returned in this buffer + + + int rtdb_next(const int handle, const int namelen, char *name) + + Return the name of the next (user inserted) entry in the data base. + The order is effectively random. + + handle = handle to RTDB + namelen = size of user provided buffer name + name = name of entry is returned in this buffer + + + int rtdb_print(const int handle, const int print_values) + + Print the contents of the data base to stdout + + handle = handle to RTDB + print_values = boolean flag ... if true values as well as + keys are printed out. + + + int rtdb_delete(const int handle, const char *name) + + Delete the entry from the database. + Return + 1 if key was present and successfully deleted + + 0 if key was not present, or if an error occured + + handle = handle to RTDB + name = name of entry to delete + +*/ + +#include + + bool rtdb_open(const char *, const char *, FILE *); + bool rtdb_close(FILE *, const char *); + bool rtdb_put(const int, const char *, const int, const int, const void *); + bool rtdb_get(const int, const char *, const int, const int, void *); + bool rtdb_get_info(const int, const char *, int *, int *, char [26]); + bool rtdb_ma_get(const int, const char *, int *, int *, int *); + bool rtdb_first(const int, const int, char *); + bool rtdb_next(const int, const int, char *); + bool rtdb_print(const int, const int); + bool rtdb_delete(const int, const char *); + +/* + Following are 'parallel' versions of the above where only + process 0 actually accesses the data base and all others + just get its output. +*/ + + bool rtdb_par_open(const char *, const char *, FILE *); + bool rtdb_par_close(FILE *, const char *); + bool rtdb_par_put(const int, const char *, const int, const int, const void *); + bool rtdb_par_get(const int, const char *, const int, const int, void *); + bool rtdb_par_get_info(const int, const char *, int *, int *, char [26]); + bool rtdb_par_ma_get(const int, const char *, int *, int *, int *); + bool rtdb_par_first(const int, const int, char *); + bool rtdb_par_next(const int, const int, char *); + bool rtdb_par_print(const int, const int); + bool rtdb_par_delete(const int, const char *); + +#endif // _RTDB_H