Began work on 'gradients' library

This commit is contained in:
Adam Parler 2023-12-06 02:53:51 -08:00 committed by Adam Parler
parent 8663962fcf
commit efd354baf1
24 changed files with 1340 additions and 119 deletions

40
Makefile Normal file
View file

@ -0,0 +1,40 @@
CC = gcc
CXX = g++
MPICC = mpicc
MPICXX = mpicxx
NWCHEM_TOP = $(shell pwd)
SRC = $(NWCHEM_TOP)/src
BIN = $(NWCHEM_TOP)/bin
BUILD = $(NWCHEM_TOP)/build
LIB_DEFINES = -DCOMPILATION_DATE="'`date +%a_%b_%d_%H:%M:%S_%Y`'" \
-DCOMPILATION_DIR="'$(TOPDIR)'" \
-DNWCHEM_BRANCH="'$(CODE_BRANCH)'"
CFLAGS=-c -Wall
LDFLAGS=
export
TARGETS=nwchem
#TARGETS := $(addprefix $(BIN)/, $(TARGETS))
.PHONY: all clean
all: $(TARGETS)
nwchem:
$(MAKE) -C $(SRC)
clean:
rm $(BIN)/$(TARGETS)
rm $(BUILD)/*.o
dist-clean: clean
rmdir $(BIN)
rmdir $(BUILD)

16
src/Makefile Normal file
View file

@ -0,0 +1,16 @@
SOURCES=
LIBRARIES=
libs: $(LIBRARY_PATH)
@mkdir -p $(LIB)
$(BIN)/nwchem: $(BUILD)/nwchem.o libs
@mkdir -p $(@D)
$(MPICC) $(LDFLAGS) $^ -o $@
$(BUILD)/%.o: $(SRC)/%.c
@mkdir -p $(@D)
$(MPICC) $(CFLAGS) -I$(SRC) -c $< -o $@

43
src/gradients/Makefile Normal file
View file

@ -0,0 +1,43 @@
# OBJ = gradients.o grad_force.o grad1.o scf_gradient.o \
grad_dens.o grad_inp.o ga_reorder.o
# OBJ_OPTIMIZE = grad2.o grad_getdens.o
# USES_BLAS = grad2.F ga_reorder.F grad_dens.F
# LIBRARY = libgradients.a
#include ../config/makefile.h
#include ../config/makelib.h
CC=gcc
BUILD = /people/parl703/nwchem/build
SRC = $(shell pwd)
LIB = /people/parl703/nwchem/lib
SOURCES := $(wildcard *.c)
# OBJECTS := $(patsubst %.c, ../../build/%.o, $(SOURCES))
# OBJ_BUILD = $(addprefix $(BUILD)/, $(OBJ) $(OBJ_OPTIMIZE))
OBJ := $(addprefix $(BUILD)/, $(OBJ))
OBJ_OPTIMIZE := $(addprefix $(BUILD)/, $(OBJ_OPTIMIZE))
all: libgradients
libgradients: object object_opt
ar -cvrsu $(LIB)/libgradients.a $(OBJ) $(OBJ_OPTIMIZE)
object: $(OBJ)
object_opt: $(OBJ_OPTIMIZE)
$(BUILD)/%.o: %.c
$(CC) -I$(SRC) -c $< -o $@

168
src/gradients/ga_reorder.c Normal file
View file

@ -0,0 +1,168 @@
#include <stdbool.h>
#include "../util/errquit.h"
//#include "global.h"
#include "../util/global.h"
void ga_reorder(int g_a, bool orow, int *rmap, bool ocol, int *cmap) {
int g_d, i, j, l_v, k_v, l_vv, k_vv, ma_type, dim1, dim2, jj;
ga_inquire(g_a, ma_type, dim1, dim2);
if (!ma_alloc_get(MT_DBL, dim1, "gareo", l_v, k_v)) {
errquit("ga_reorder: could not allocate column", dim1, MA_ERR);
}
if (!ma_alloc_get(MT_DBL, dim1, "gareo", l_vv, k_vv)) {
errquit("ga_reorder: could not allocate column2", dim1, MA_ERR);
}
ga_sync();
if (!ga_duplicate(g_a, g_d, "ga_reorder")) {
errquit("ga_reorder: duplicate failed", 0, GA_ERR);
}
ga_copy(g_a, g_d);
for (j = ga_nodeid(); j < dim2; j += ga_nnodes()) {
if (orow) {
ga_get(g_d, 1, dim1, j, j, dbl_mb[k_v], dim1);
for (i = 0; i < dim1; i++) {
dbl_mb[k_vv+rmap[i]-1] = dbl_mb[k_v+i-1];
}
} else {
ga_get(g_d, 1, dim1, j, j, dbl_mb[k_vv], dim1);
}
jj = j;
if (ocol) {jj = cmap[j];}
ga_put(g_a, 1, dim1, jj, jj, dbl_mb[k_vv], dim1);
}
if (!ma_free_heap(l_vv)) {
errquit("ga_reo: ma?", 0, MA_ERR);
}
if (!ma_free_heap(l_v)) {
errquit("ga_reo: ma2?", 0, MA_ERR);
}
if (!ga_destroy(g_d)) {
errquit("ga_reo: ga_destroy?", 0, GA_ERR);
}
}
void nga_reorder(int g_a, bool orow, int *rmap, bool ocol, int *cmap) {
/*
This is basically just an extension of ga_reorder and is not very
generic at this point. As a matter of fact, it assumes (and tests)
that the dimension is 3 and that you only want to reorder the last
two indices. This can be made more general after I test this version.
Also, I am wasting a lot of memory by duplicating the whole ga. This
will need to be optimized in the future.
*/
int g_d, i, j, l_v, k_v, l_vv, k_vv, ma_type;
int dim0, dim1, dim2, jj;
int ndim, dims[3], lo[3], hi[3], ld[2];
ndim = ga_ndim(g_a);
if (ndim != 3) {
errquit("nga_reorder: must have 3 dimensions", ndim, GA_ERR);
}
nga_inquire(g_a, ndim, dims);
dim1 = dims[1];
dim2 = dims[2];
if (!ma_alloc_get(MT_DBL, dim1, "gareo", l_v, k_v)) {
errquit("ga_reorder: could not allocate column", dim1, GA_ERR);
}
if (!ma_alloc_get(MT_DBL, dim1, "gareo", l_vv, k_vv)) {
errquit("ga_reorder: could not allocate column2", dim1, GA_ERR);
}
ga_sync();
if (!ga_duplicate(g_a, g_d, "ga_reorder")) {
errquit("ga_reorder: duplicate failed", 0, GA_ERR);
}
ga_copy(g_a, g_d);
ld[0] = 1;
lo[1] = 1;
hi[1] = dim1;
ld[1] = dim1;
for (dim0 = 0; dim0 < dims[0]; i++) {
lo[0] = dim0;
hi[0] = dim0;
for (j = ga_nodeid(); j < dim2; j += ga_nnodes()) {
lo[2] = j;
hi[2] = j;
if (orow) {
nga_get(g_d, lo, hi, dbl_mb[k_v], ld);
for (i = 0; i < dim1; i++) {
dbl_mb[k_vv+rmap[i]-1] = dbl_mb[k_v+i-1];
}
} else {
nga_get(g_d, lo, hi, dbl_mb[k_vv], ld);
}
jj = j;
if (ocol) {jj = cmap[j];}
lo[2] = jj;
hi[2] = jj;
nga_put(g_a, lo, hi, dbl_mb[k_vv], ld);
}
}
if (!ma_free_heap(l_vv)) errquit("ga_reo: ma?", 0, MA_ERR);
if (!ma_free_heap(l_v)) errquit("ga_reo: ma2?", 0,MA_ERR);
ga_sync();
if (!ma_free_heap(g_d)) errquit("ga_reo: ga_destroy", 0, GA_ERR);
}
void matrix_reorder(int dim1, int dim2, double *a, bool orow, int *rmap, bool ocol, int *cmap) {
int i, j, l_v, k_v, l_vv, k_vv, jj;
int l_d, k_d;
if (!ma_alloc_get(MT_DBL, dim1*dim2, "mareo", l_d, k_d)) {
errquit("ga_reorder: could not allocate dup", dim1*dim2, MA_ERR);
}
if (!ma_alloc_get(MT_DBL, dim1, "mareo", l_v, k_v)) {
errquit("ga_reorder: could not allocate column", dim1, MA_ERR);
}
if (!ma_alloc_get(MT_DBL, dim1, "mareo", l_vv, k_vv)) {
errquit("ga_reorder: could not allocate column2", dim1, MA_ERR);
}
dcopy(dim1*dim2, a, 1, dbl_mb[k_d], 1);
for (j = 0; j < dim2; j++) {
if (orow) {
dcopy(dim1, dbl_mb[k_d+j*dim1], 1, dbl_mb[k_v], 1);
for (i = 0; i < dim1; i++) {
dbl_mb[k_vv+rmap[i]-1] = dbl_mb[k_v+i-1];
}
} else {
dcopy(dim1, dbl_mb[k_d+j*dim1], 1, dbl_mb[k_vv], 1);
}
jj = j;
if (ocol) jj = cmap[j];
dcopy(dim1, dbl_mb[k_vv], 1, a[jj], 1);
}
if (!ma_free_heap(l_vv)) {
errquit("ma_reo: ma?", 0, MA_ERR);
}
if (!ma_free_heap(l_v)) {
errquit("ma_reo: ma2?", 0, MA_ERR);
}
if (!ma_free_heap(l_d)) {
errquit("ma_reo: ma?", 0, MA_ERR);
}
}

218
src/gradients/grad1.c Normal file
View file

@ -0,0 +1,218 @@
#include <stdbool.h>
#include "../util/global.h"
//#include "geom.h"
//#include "bas.f"
//#include "rtdb.h"
//#include "sym.h"
//#include "bq_params.h"
#define NO_BQGEM 1
void grad1(double *H, int lbuf, double *scr, int lscr, double *dens,
double *wdens, double *frc_nuc, double *frc_kin, double *frc_wgh,
int g_force, int *g_dens, int g_wdens, int basis, int geom, int nproc,
int nat, int max_at_bf, int rtdb, bool oskel, int ndens ) {
int ijatom, next, iat1, iat2, iat3, ish1, ish2,
iab1f, iab1l, iab2f, iab2l, iac1f, iac1l, iac2f, iac2l,
if1, il1, if2, il2, icart, ic, nint, ip1, ip2;
double crd1[3], crd2[3]; // atomic coordinates;
int idatom[2];
double dE, dx, dy, dz, qfac, fact, q1, q2;
bool status, pointforce, dobq;
char name[16];
int bq_ncent;
int i_qbq,i_cbq;
double r12;
int task_size;
// AJL/Begin/SPIN ECPs
int ecp_channels;
int iecp;
double H_beta[lbuf];
double dens_beta[max_at_bf][max_at_bf];
#ifdef NO_BQGEM
//#include "inp.h"
char bqchar[2];
#endif
// Read this value from rtdb vvvv
if (!rtdb_get(rtdb, "dft:spin_polarised_ecps'", MT_INT, 1, ecp_channels)) {
ecp_channels = 1;
}
/* AJL: With spin-polarised ECPs Hcore will be spin dependent
See Szabo and Ostlund pg. 215
So we need to separate out the densities
if (ecp_channels.gt.1) then
Restore alpha and beta densities to calculate spin-polarised
derivatives
call ga_print(g_dens(1))
call ga_print(g_dens(2))
call ga_dadd(1d0, g_dens(1), -1d0, g_dens(2), g_dens(1))
call ga_print(g_dens(1))
call ga_print(g_dens(2))
end if
AJL/End */
task_size = 1;
status = rtdb_parallel(true); // Broadcast reads to all processes
pointforce = geom_include_bqbq(geom);
dobq = geom_extbq_on();
hf_print_set(1);
ijatom = -1;
next = nxtask(nproc,task_size);
for (iat1 = 0; iat1 < nat; iat1++) {
for (iat2 = 0; iat2 < iat1; iat2++) {
ijatom++;
if (ijatom == next) {
status = bas_ce2bfr(basis,iat1,iab1f,iab1l);
status = bas_ce2bfr(basis,iat2,iab2f,iab2l);
if (iab1f <= 0 || iab2f <= 0) {
// At least one center has no functions on it ... next atom
goto g1010;
}
if (oskel) {
if (!sym_atom_pair(geom, iat1, iat2, qfac)) goto g1010;
} else {
qfac = 1.0;
}
status = bas_ce2cnr(basis,iat1,iac1f,iac1l);
status = bas_ce2cnr(basis,iat2,iac2f,iac2l);
// AJL/Begin/SPIN ECPs
// call ga_get(g_dens,iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
for (iecp = 0; iecp < ecp_channels; iecp++) {
if (iecp == 1) {
ga_get(g_dens[iecp],iab1f,iab1l, iab2f,iab2l,dens,max_at_bf);
} else {
ga_get(g_dens[iecp],iab1f,iab1l, iab2f,iab2l,dens_beta,max_at_bf);
}
}
// Recombine g_dens, as it is not used again
// if (ecp_channels.gt.1) then
// call ga_dadd(1d0, g_dens(1), 1d0, g_dens(2), g_dens(1))
// end if
// g_wdens is not dependent on spin, so can leave this
ga_get(g_wdens,iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf);
// AJL/End
for (ish1 = iac1f; ish1 < iac1l; ish1++) {
if ( iat1 == iat2 ) iac2l = ish1;
for (ish2 = iac2f; ish2 < iac2l; ish2++) {
// shell block in atomic (D/Dw)-matrix block
status = bas_cn2bfr(basis,ish1,if1,il1);
if1 = if1 - iab1f + 1;
il1 = il1 - iab1f + 1;
status = bas_cn2bfr(basis,ish2,if2,il2);
if2 = if2 - iab2f + 1;
il2 = il2 - iab2f + 1;
nint = ( il1 - if1 + 1 ) * ( il2 - if2 + 1 );
// overlap derivatives
intd_1eov(basis,ish1,basis,ish2,lscr,scr, lbuf,H,idatom);
// Dw x S
if ( idatom[0] >= 1 ) {
// idatom(1).ge.0 <=> idatom(2).ge.0 (no check necessary)
ic = 0;
for (icart = 0; icart < 3; icart++) {
dE = 0.0;
for (ip1 = if1; ip1 < il1; ip1++) {
for (ip2 = if2; ip2 < il2; ip2++) {
dE += wdens[ip1*il1+ip2] * H[ic];
}
}
dE = dE * qfac;
frc_wgh[3*icart+idatom[0]] = frc_wgh[3*icart+idatom[0]] - dE - dE;
frc_wgh[3*icart+idatom[1]] = frc_wgh[3*icart+idatom[1]] + dE + dE;
}
}
// 1el. derivatives
if (!dobq) {
intd_1eh1(basis,ish1,basis,ish2,lscr,scr,lbuf,H);
} else {
intd_1epot(basis,ish1,basis,ish2,lscr,scr,lbuf,H);
}
// AJL/Begin/SPIN ECPs
// With spin-polarised ECPs Hcore will be spin dependent
// See Szabo and Ostlund pg. 215
if (ecp_channels > 1) {
// 1el. derivatives
if (!dobq) {
// For now this will do, but this could be more efficiently done
intd_1eh1_beta(basis,ish1,basis,ish2,lscr,scr,lbuf,H_beta);
} else {
intd_1epot_beta(basis,ish1,basis,ish2,lscr,scr,lbuf,H_beta);
}
}
// AJL/End
// D x H
ic = 0;
for (iat3 = 0; iat3 < nat; iat3++) {
for (icart = 0; icart < 3; icart++) {
dE = 0.0;
}
}
}
}
}
g1010: continue;
}
}
}
/*
C> \brief calculate the gradient terms due to the interaction with the
C> COSMO charges
C>
C> Evaluate the gradient contributions from the COSMO embedding. The
C> original part is from Klamt and Sch&uuml;&uuml;rmann [1]
C> (see Eqs.(13-16)). The derivatives of matrix \f$A\f$ have been
C> modified by York and Karplus [2] (see Eqs.(73-76)) to obtain smooth
C> potential energy surfaces. York and Karplus also modified matrix
C> \f$B\f$ which is easy to do in their classical force field code.
C> In an ab-initio code this not so easy to do and as it is not
C> required to eliminate singularities the original expression from [1]
C> for \f$B\f$ is used here.
C>
C> ### References ###
C>
C> [1] A. Klamt, G. Sch&uuml;&uuml;rmann,
C> "COSMO: a new approach to dielectric screening in solvents with
C> explicit expressions for the screening energy and its gradient",
C> <i>J. Chem. Soc., Perkin Trans. 2</i>, 1993, pp 799-805, DOI:
C> <a href="https://doi.org/10.1039/P29930000799">
C> 10.1039/P29930000799</a>.
C>
C> [2] D.M. York, M. Karplus,
C> "A smooth solvation potential based on the conductor-like
C> screening model", <i>J. Phys. Chem. A</i> (1999) <b>103</b>,
C> pp 11060-11079, DOI:
C> <a href="https://doi.org/10.1021/jp992097l">
C> 10.1021/jp992097l</a>.
*/

0
src/gradients/grad2.c Normal file
View file

View file

View file

View file

0
src/gradients/grad_inp.c Normal file
View file

View file

128
src/gradients/gradients.c Normal file
View file

@ -0,0 +1,128 @@
#include <stdio.h>
#include <string.h>
#include <stdbool.h>
#include "../util/errquit.h"
//#include "bas.h"
//#include "geom.h"
#include "../util/global.h"
#include "../rtdb/rtdb.h"
//#include "schwarz.h"
#include "../util/util.h"
#include "../util/stdio.h"
bool gradients(int rtdb) {
int mt_log;
// gradients module.
/*
Assumes SCF has been completed, MO vectors stored
and all information is still in the RTDB
*/
int geom, basis; // handles
bool status;
char title[255];
bool odbug;
bool ocosmo;
bool osome;
status = rtdb_parallel(true); // Broadcast reads to all processes
ecce_print_module_entry("gradients");
// Extract high level info from the data-base setting defaults
if (!rtdb_cget(rtdb, "title", 1, title)) strncpy(title, " ", 4);
if (!geom_create(geom, "geometry")) errquit("gradients: geom_create?", 0, GEOM_ERR);
if (!geom_rtdb_load(rtdb, geom, "geometry")) errquit("gradients: no geometry ", 0, GEOM_ERR);
if (!bas_create(basis, "ao basis")) errquit("gradients: bas_create?", 0, BASIS_ERR);
if (!bas_rtdb_load(rtdb, geom, basis, "ao basis")) errquit("gradients: no ao basis", 0, BASIS_ERR);
if (!int_normalize(rtdb,basis)) errquit("gradients: normalization failed", 911, INT_ERR);
/*
Figure out the numer of electrons from the required total
charge and the sum of nuclear charges
if (.not. rtdb_get(rtdb, 'charge', MT_DBL, 1, charge))
$ charge = 0.0d0
*/
if (nodeid == 0) {
if (util_print("information", print_low)) {
util_print_centered(LuOut, "NWChem Gradients Module", 40, true);
fprintf(stdout, "%s", LuOut);
util_flush();
}
if (util_print("information", print_medium)) {
fprintf(stdout, "%s", LuOut);
if (title != " ") {
util_print_centered(LuOut, title, 40, false);
fprintf(stdout, "%s", LuOut);
}
util_flush(LuOut);
}
if (util_print("geometry", print_high)) {
if (!geom_print(geom)) {
errquit("gradients: geom_print ?", 0, GEOM_ERR);
}
util_flush(LuOut);
}
if (uitl_print("basis", print_high)) {
if (!bas_print(basis)) {
errquit("gradients: bas_print ?", 0, BASIS_ERR);
}
util_flush(LuOut);
}
}
odbug = false;
odbug = odbug && ga_nodeid() == 0;
if (rtdb_get(rtdb,"slv:cosmo", mt_log, 1, ocosmo)) {
if (odbug) {
fprintf(stdout, "-cosmo- ... found in -gradients-%s %d",
ocosmo ? "true" : "false", ga_nodeid());
}
if (ocosmo) {
if (odbug) {
osome = true;
} else {
osome = false;
}
osome = osome && ga_nodeid() == 0;
if (odbug) {
fprintf(stdout, "-cosmo- ... found and .true. %s %d",
ocosmo ? "true" : "false", ga_nodeid());
}
} else {
if (odbug) {
fprintf(stdout, "-cosmo- ... found but .false. %s %d",
ocosmo ? "true" : "false", ga_nodeid());
}
}
} else {
if (odbug) {
fprintf(stdout, "-cosmo- not found in -gradients-");
}
}
ga_sync();
// go for it ... finally ...
grad_force(rtdb, basis, geom);
// gradients is done destroy basis and geometry handles
// (e.g., preserve the memory available to other modules!!)
if ( !(bas_destroy(basis) && geom_destroy(geom)) ) {
errquit("gradients:error destroying geom and basis handles",911, GEOM_ERR);
}
ecce_print_module_exit("gradients","ok");
return true;
}

View file

@ -0,0 +1,35 @@
#include <stdbool.h>
#include "../rtdb/rtdb.h"
#include "../util/errquit.h"
bool mcscf_gradient(int rtdb) {
if (!mcscf(rtdb)) {
errquit("mcscf_gradient: mcscf energy failed", 0, CALC_ERR);
}
util_print_push();
util_print_rtdb_load(rtdb,"mcscf");
if(!gradients(rtdb)) {
errquit("mcscf_gradient: gradients failed", 0, CALC_ERR);
}
util_print_pop();
return true;
}
bool scf_gradient(int rtdb) {
if (!scf(rtdb)) {
errquit("scf_gradient: scf energy failed", 0, CALC_ERR);
}
util_print_push();
util_print_rtdb_load(rtdb, "scf");
if (!gradients(rtdb)) {
errquit("scf_gradient: gradients failed", 0, CALC_ERR);
}
util_print_pop();
return true;
}

View file

@ -28,50 +28,51 @@
int32_t IO_CODE;
#endif
// $Id$
/*
$Id$
// ======================================================================================================
//> \mainpage Northwest Computational Chemistry Package (NWChem) 7.0.2
//>
//> NWChem is an open-source computational chemistry package distributed under the terms of
//> the Educational Community License (ECL) 2.0
//>
//> This software and its documentation were developed at the EMSL at Pacific Northwest National Laboratory,
//> a multiprogram national laboratory, operated for the U.S. Department of Energy by Battelle under
//> Contract Number DE-AC05-76RL01830. Support for this work was provided by the Department of Energy
//> Office of Biological and Environmental Research, Office of Basic Energy Science, and the Office of
//> Advanced Scientific Computing.
//>
//> Licensed under the Educational Community License, Version 2.0 (the "License"); you may
//> not use this file except in compliance with the License. You may obtain a copy of the
//> License at <a href="https://opensource.org/licenses/ECL-2.0">https://opensource.org/licenses/ECL-2.0</a>.
//>
//> Unless required by applicable law or agreed to in writing, software distributed under the
//> License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
//> either express or implied. See the License for the specific language governing
//> permissions and limitations under the License.
//>
//> Further information, including user documentation and forums, may be found at
//> <a href="http://www.nwchem-sw.org/">http://www.nwchem-sw.org/</a>. Alternatively,
//> the paper
//>
//> * M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, T.P. Straatsma, H.J.J. Van Dam,
//> D. Wang, J. Nieplocha, E. Apra, T.L. Windus, W.A. de Jong (2010)<br>
//> "NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations"<br>
//> <I>Computer Physics Communications</I>, <b>181</b>, 14771489, DOI: <a href="https://doi.org/10.1016/j.cpc.2010.04.018">10.1016/j.cpc.2010.04.018</a>
//>
//> provides details on the codes capabilities.
//>
//> Copyright (c) 1994-2020 Pacific Northwest National Laboratory, Battelle Memorial Institute
//>
//> Environmental Molecular Sciences Laboratory (EMSL)<br>
//> Pacific Northwest National Laboratory<br>
//> Richland, WA 99352
=====================================================================================================
\mainpage Northwest Computational Chemistry Package (NWChem) 7.0.2
// ======================================================================================================
NWChem is an open-source computational chemistry package distributed under the terms of
the Educational Community License (ECL) 2.0
This software and its documentation were developed at the EMSL at Pacific Northwest National Laboratory,
a multiprogram national laboratory, operated for the U.S. Department of Energy by Battelle under
Contract Number DE-AC05-76RL01830. Support for this work was provided by the Department of Energy
Office of Biological and Environmental Research, Office of Basic Energy Science, and the Office of
Advanced Scientific Computing.
int main(){
Licensed under the Educational Community License, Version 2.0 (the "License"); you may
not use this file except in compliance with the License. You may obtain a copy of the
License at <a href="https://opensource.org/licenses/ECL-2.0">https://opensource.org/licenses/ECL-2.0</a>.
Unless required by applicable law or agreed to in writing, software distributed under the
License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
either express or implied. See the License for the specific language governing
permissions and limitations under the License.
Further information, including user documentation and forums, may be found at
<a href="http://www.nwchem-sw.org/">http://www.nwchem-sw.org/</a>. Alternatively,
the paper
* M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, T.P. Straatsma, H.J.J. Van Dam,
D. Wang, J. Nieplocha, E. Apra, T.L. Windus, W.A. de Jong (2010)<br>
"NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations"<br>
<I>Computer Physics Communications</I>, <b>181</b>, 14771489, DOI: <a href="https://doi.org/10.1016/j.cpc.2010.04.018">10.1016/j.cpc.2010.04.018</a>
provides details on the codes capabilities.
Copyright (c) 1994-2020 Pacific Northwest National Laboratory, Battelle Memorial Institute
Environmental Molecular Sciences Laboratory (EMSL)<br>
Pacific Northwest National Laboratory<br>
Richland, WA 99352
=====================================================================================================
*/
int main(int argc, char *argv[]) {
char input_filename[nw_max_path_len], rtdb_name[nw_max_path_len];
double total_wall, total_cpu;
#ifdef USE_OFFLOAD
@ -150,20 +151,21 @@ int main(){
errquit('nwchem.F: ma_init failed (ga_uses_ma=F)',911, MA_ERR);
}
}
/*
Touch OpenMP here so that any runtime initialization happens up-front.
This ensures that any printout that the OpenMP runtime generates,
such as affinity information, appears at the top of the output file.
Otherwise, it might not appear until e.g. the CCSD module, at which
point it will pollute the output file in an undesirable way.
// Touch OpenMP here so that any runtime initialization happens up-front.
// This ensures that any printout that the OpenMP runtime generates,
// such as affinity information, appears at the top of the output file.
// Otherwise, it might not appear until e.g. the CCSD module, at which
// point it will pollute the output file in an undesirable way.
Do not move this in front of GA/MPI/TCGMSG initialization, since the
OpenMP runtime may inherit affinity information from MPI that is only
determined during MPI initialization.
// Do not move this in front of GA/MPI/TCGMSG initialization, since the
// OpenMP runtime may inherit affinity information from MPI that is only
// determined during MPI initialization.
// Format definition is outside of preprocessor protection to ensure the
// label is not accidentally reused, since that will not be caught by
// testing that does not enable OpenMP.
Format definition is outside of preprocessor protection to ensure the
label is not accidentally reused, since that will not be caught by
testing that does not enable OpenMP.
*/
g99 format(2x, 'NWChem w/ OpenMP: Maximum threads = ',i4);
#if defined(USE_OPENMP){
@ -192,27 +194,27 @@ int main(){
errquit('nwchem.cpp: ma_set_numalign failed', 911, MA_ERR);
}
//old:------------------------------------------------------- START ---------
//old:C GA allocations come out of MA space, so lump them together
//old:C and let MA impose the limits on GA sizes instead of actually
//old:C using the global limit.
//old:C
//old: if ( ga_uses_ma() ) then
//old: if (.not. ma_init(mt_dbl, stack, heap+global))
//old: $ call errquit('nwchem: ma_init failed', -1)
//old: call ga_initialize
//old:C
//old:C GA allocations are separate from MA, so the separate limit
//old:C must be enforced. Note GA only understands bytes.
//old:C
//old: else
//old: if (.not. ma_init(mt_dbl, stack, heap))
//old: $ call errquit('nwchem: ma_init failed', -1)
//old: call ga_initialize_ltd(ma_sizeof(mt_dbl, global, mt_byte) )
//old: endif
//old:------------------------------------------------------- END -----------
/*
old:------------------------------------------------------- START ---------
old:C GA allocations come out of MA space, so lump them together
old:C and let MA impose the limits on GA sizes instead of actually
old:C using the global limit.
old:C
old: if ( ga_uses_ma() ) then
old: if (.not. ma_init(mt_dbl, stack, heap+global))
old: $ call errquit('nwchem: ma_init failed', -1)
old: call ga_initialize
old:C
old:C GA allocations are separate from MA, so the separate limit
old:C must be enforced. Note GA only understands bytes.
old:C
old: else
old: if (.not. ma_init(mt_dbl, stack, heap))
old: $ call errquit('nwchem: ma_init failed', -1)
old: call ga_initialize_ltd(ma_sizeof(mt_dbl, global, mt_byte) )
old: endif
old:------------------------------------------------------- END -----------
*/
//*** call nxtval_ga_initialize()
// Trap SIGFPE after GA to override handler
@ -401,3 +403,148 @@ int main(){
return 0;
}
/*
void nwchem_head_info(char* argv[]) {
int ierr, num_procs, nodeid;
ierr = MPI_Init(&argc, &argv);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &nodeid);
ierr = MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
ierr = MPI_Finalize();
FILE *fptr;
time_t timer;
struct tm* tm_info;
char compiled[] = __TIMESTAMP__;
char nwchem_rev[] = VERSION;
char branch[] = NWCHEM_BRANCH;
char *raw_srcdir = realpath(argv[0], NULL);
char c, hostname[80], executable[nw_max_path_len], date[26], input_filename[nw_max_path_len];
char ga_rev[nw_max_path_len], srcdir[nw_max_path_len], thing[32][nw_max_path_len], *ptr;
char rtdb_name[nw_max_path_len], file_prefix[nw_max_path_len], folder_prefix[nw_max_path_len];
char cstart[10];
int host_return, i, depth, nproc;
strncpy(input_filename, argv[1], nw_max_path_len-1);
printf(" argument 1 = %s\n\n", input_filename);
// printf("%.*s", 30, "=================");
//printf("%0*d\n", 20, 0);
printf("\n\n==============================");
printf(" echo of input deck ");
printf("==============================\n");
// Open file
fptr = fopen(input_filename, "r");
if (fptr == NULL)
{
printf("Cannot open file \n");
exit(0);
}
// Read contents from file
c = fgetc(fptr);
while (c != EOF)
{
printf("%c", c);
c = fgetc(fptr);
}
// Close file
fclose(fptr);
printf("\n==========================================");
printf("======================================\n\n\n\n\n\n\n");
// Printing hostname
host_return = gethostname(hostname, sizeof(hostname));
if (host_return == -1) errquit("nwchem: failed to get hostname", 0, 10);
// Printing program
strncpy(executable, argv[0], nw_max_path_len-1);
// Printing current date
timer = time(NULL);
tm_info = localtime(&timer);
strftime(date, 30, "%a %b %d %H:%M:%S %Y", tm_info);
// Reformatting compilation date
ptr = compiled;
while (*ptr) {
if (*ptr == ' ')
*ptr = '_';
ptr++;
}
// Getting top-level source folder
depth = 0;
ptr = strtok(raw_srcdir, "/");
while (ptr != NULL) {
strcpy(thing[depth], ptr);
ptr = strtok(NULL, "/");
depth++;
}
depth -= 3;
i = 0;
srcdir[0] = '\0';
while (i < depth) {
strcat(srcdir, "/");
strcat(srcdir, thing[i+1]);
i++;
}
// Getting release info (OLD)
/*#ifdef RELEASE
#define NWCHEM_BRANCH "7.0.2"
#else
#define NWCHEM_BRANCH "Development"
#endif/*
// Printing GA info
strncpy(ga_rev, "5.7.2", 79);
snprintf(file_prefix, 79, "%s.", "eu_hdehp_cmpx");
strncpy(folder_prefix, "./perm", nw_max_path_len-2);
snprintf(rtdb_name, nw_max_path_len, "%s/%sdb", folder_prefix, file_prefix);
strncpy(cstart, "startup", 9);
#if defined(MPI)
MPI_Comm_size(MPI_COMM_WORLD, &nproc);
#elif defined(_OPENMP)
nproc = omp_get_num_threads();
#else
nproc = 1;
#endif
printf(" Job information\n");
printf(" ---------------\n");
printf(" hostname = %s\n", hostname);
printf(" program = %s\n", executable);
printf(" date = %s\n\n", date);
printf(" compiled = %s\n", compiled);
printf(" source = %s\n", srcdir);
printf(" nwchem branch = %s\n", branch);
printf(" nwchem revision = %s\n", nwchem_rev);
printf(" ga revision = %s\n", ga_rev);
printf(" use scalapack = %s\n", util_scalapack_info() ? "T" : "F");
printf(" input = %s\n", input_filename);
printf(" prefix = %s\n", file_prefix);
printf(" data base = %s\n", rtdb_name);
printf(" status = %s\n", cstart);
printf(" nproc = %8d\n", nproc);
printf(" time left = %6ds\n", util_batch_job_time_remaining());
*/

39
src/rtdb/Makefile Normal file
View file

@ -0,0 +1,39 @@
LIBRARY = libnwcutil.a
LIBRARIES += $(LIBRARY)
OBJ_OPTIMIZE += rtdb.o rtdb_seq.o context.o
HEADERS = context.h rtdb.h rtdb.cray.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 testgr.o testgr
TEST_LIBS = $(LIBRARY) $(LIBS)
$(LIBRARY): $(LIB_TARGETS)
davetest: davetest.o $(LIBRARY_PATH)
$(CC) $(CFLAGS) $(LDFLAGS) -o $@ davetest.o $(LIBS)
cntx: cntx.o $(LIBRARY_PATH)
$(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ $(LIBS)
interact: interact.o $(LIBRARY_PATH)
$(CC) $(CFLAGS) $(LDFLAGS) -o $@ interact.o $(LIBRARY_PATH) -lglobal -ltcgmsg -lm
test: test.o $(LIBRARY_PATH)
$(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ $(LIBS)
testgr: testgr.o $(LIBRARY_PATH)
$(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ $(LIBS)
rtdbtest: rtdbtest.o $(LIBRARY_PATH)
$(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ -lm $(LIBS)
rtdbpartest: rtdb_par_test.o $(LIBRARY)
$(CC) $(CFLAGS) -o $@ rtdb_par_test.o $(TEST_LIBS)
context: context.o $(LIBRARY)
$(CC) $(CFLAGS) -o $@ context.o $(TEST_LIBS)

210
src/rtdb/context.c Normal file
View file

@ -0,0 +1,210 @@
/*$Id$*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "rtdb.h"
#include "macdecls.h"
#include "misc.h"
#define MAX_CLEN 4096
static char context[MAX_CLEN];
int context_set(const char *string)
{
if (strlen(string) < sizeof(context)) {
(void) strcpy(context, string);
return 1;
}
else {
fprintf(stderr, "context_set: string too long? %s\n", string);
fflush(stderr);
return 0;
}
}
char *context_get(void)
{
return strdup(context);
}
int context_rtdb_store(int rtdb)
{
return rtdb_put(rtdb, "Context", MT_CHAR, strlen(context)+1, context);
}
int context_rtdb_load(int rtdb)
{
return rtdb_get(rtdb, "Context", MT_CHAR, sizeof(context), context);
}
int context_push(const char *string)
{
int clen = strlen(context);
int slen = strlen(string);
if (slen+clen+2 >= sizeof(context)) {
fprintf(stderr, "context_push: static dimension of context too small\n");
fprintf(stderr, "context_push: current = %s\n", context);
fprintf(stderr, "context_push: pushing = %s\n", string);
return 0;
}
else {
(void) strcpy(context+clen, string);
(void) strcpy(context+clen+slen, ":");
return 1;
}
}
int context_pop(const char *string)
{
int clen = strlen(context);
int slen = strlen(string);
if (clen)
clen--; /* Trailing colon */
if (slen <= clen && strncmp(context+clen-slen, string, slen) == 0) {
context[clen-slen] = 0;
return 1;
}
else {
fprintf(stderr, "context_pop: current = %s\n", context);
fprintf(stderr, "context_pop: popping = %s\n", string);
return 0;
}
}
int context_rtdb_match(int rtdb, const char *name, int reslen,
char *result)
{
char buf[MAX_CLEN];
int blen = strlen(context);
if (blen+strlen(name)+1 > sizeof(buf)) {
fprintf(stderr, "context_rtdb_match: buffer size exceeded\n");
fprintf(stderr, "context_rtdb_match: current = %s\n", context);
fprintf(stderr, "context_rtdb_match: pushing = %s\n", name);
return 0;
}
strcpy(buf, context);
while (1) {
int ma_type, nelem;
char date[26];
/* Append name to current context */
(void) strcpy(buf+blen, name);
if (rtdb_get_info(rtdb, buf, &ma_type, &nelem, date)) {
if (ma_type == MT_CHAR) {
if (!rtdb_get(rtdb, buf, ma_type, reslen, result)) {
fprintf(stderr, "context_rtdb_match: rtdb_get failed?\n");
return 0;
}
reslen = strlen(result);
if (result[reslen-1] == '\n') /* Fortran cput appends an unwanted CR */
result[reslen-1] = 0;
return 1;
}
else {
fprintf(stderr, "context_rtdb_match: found %s but is wrong type\n",
name);
return 0;
}
}
else {
/* Did not find entry ... pop the context stack */
if (!blen)
return 0; /* Stack is alredy empty */
blen--;
while (--blen > 0)
if (buf[blen] == ':')
break;
}
}
return 1; /* Never executed */
}
int context_prefix(const char *name, char *result, int result_len)
{
if ((strlen(name)+strlen(context)+1) > result_len) {
fprintf(stderr, "constant_prefix: result too short\n");
return 0;
}
strcpy(result,context);
strcpy(result+strlen(context),name);
return 1;
}
/*
static void context_print()
{
printf("context = -%s-\n", context);
}
int main()
{
int rtdb;
char *cntx;
(void) MA_initialize(MT_CHAR, -1, -1);
if (!rtdb_open("test.db", "unknown", &rtdb))
error("testcontext: open failed on %s\n", "test.db");
context_print();
if (!context_push("optimize"))
error("context push failed %d\n", 0);
context_print();
if (!context_push("scf"))
error("context push failed %d\n", 0);
context_print();
if (!context_push("rhf"))
error("context push failed %d\n", 0);
context_print();
if (!context_push("pcg"))
error("context push failed %d\n", 0);
context_print();
(void) context_store(rtdb);
(void) context_set("");
(void) context_print();
if (!context_load(rtdb))
error("context_load: failed %d\n", 0);
(void) context_print();
cntx = context_get();
printf("context from get = %s\n", cntx);
if (context_pop("scf"))
error("context pop succeeded %d\n", 0);
if (!context_pop("pcg"))
error("context pop failed %d\n", 0);
context_print();
if (!context_pop("rhf"))
error("context pop failed %d\n", 0);
context_print();
if (!context_pop("scf"))
error("context pop failed %d\n", 0);
context_print();
if (!context_pop("optimize"))
error("context pop failed %d\n", 0);
context_print();
(void) rtdb_close(rtdb, "delete");
return 0;
}
*/

14
src/rtdb/context.h Normal file
View file

@ -0,0 +1,14 @@
/*$Id$*/
int context_set(const char *);
char *context_get(void);
int context_rtdb_store(int);
int context_rtdb_load(int);
int context_push(const char *);
int context_pop(const char *);
int context_rtdb_match(int, const char *, int, char *);
int context_prefix(const char *, char *, int);
#if defined(CRAY) || defined(WIN32)
#include "rtdb.cray.h"
#endif

15
src/rtdb/davetest.c Normal file
View file

@ -0,0 +1,15 @@
int main(int argc, char *argv[]) {
char name[128];
int rtdb;
int crap;
strncpy(name, "h2o.db", 12);
pbeinf();
return 0;
}

26
src/rtdb/rtdb.cray.h Normal file
View file

@ -0,0 +1,26 @@
/*$Id$*/
#if (defined(CRAY) || defined(WIN32)) &&!defined(__crayx1) &&!defined(__MINGW32__)
#define context_pop_ CONTEXT_POP
#define context_prefix_ CONTEXT_PREFIX
#define context_push_ CONTEXT_PUSH
#define context_rtdb_load_ CONTEXT_RTDB_LOAD
#define context_rtdb_match_ CONTEXT_RTDB_MATCH
#define context_rtdb_store_ CONTEXT_RTDB_STORE
#define context_set_ CONTEXT_SET
#define context_get_ CONTEXT_GET
#define rtdb_cget_ RTDB_CGET
#define rtdb_close_ RTDB_CLOSE
#define rtdb_cput_ RTDB_CPUT
#define rtdb_delete_ RTDB_DELETE
#define rtdb_first_ RTDB_FIRST
#define rtdb_get_ RTDB_GET
#define rtdb_get_info_ RTDB_GET_INFO
#define rtdb_ma_get_ RTDB_MA_GET
#define rtdb_next_ RTDB_NEXT
#define rtdb_open_ RTDB_OPEN
#define rtdb_parallel_ RTDB_PARALLEL
#define rtdb_put_ RTDB_PUT
#define rtdb_print_ RTDB_PRINT
#define rtdb_print_usage_ RTDB_PRINT_USAGE
#endif

32
src/rtdb/testgr.c Normal file
View file

@ -0,0 +1,32 @@
#include <stdbool.h>
#include "rtdb.h"
int main(int argc, char *argv[]) {
int rtdb, ma_handle, ma_index;
int itest[3], ibuf[3];
float ftest[4], fbuf[4];
double dtest[5], dbuf[5];
char cbuf[4][20], ccbuf[4][20];
char name[20], rtdb_fname[20];
char date[26];
bool status;
int type, nelem, i;
itest = {1, 2, 3};
ftest = {1.0, 2.0, 3.0, 4.0};
dtest = {1.0, 2.0, 3.0, 4.0, 5.0};
cbuf[0] = "Have";
cbuf[1] = "a";
cbuf[2] = "nice";
cbuf[3] = "day, Robert!";
pbeginf();
if (!ma_init(MT_DBL, -1, -1)) exit;
ga_initialize()
}

11
src/util/global.h Normal file
View file

@ -0,0 +1,11 @@
#ifndef _GLOBAL_H
#define _GLOBAL_H
const int MT_DBL = 8;
const int MT_INT = 4;
int nodeid();
double *dbl_mb;
#endif

12
src/util/printlevels.h Normal file
View file

@ -0,0 +1,12 @@
#ifndef _PRINTLEVELS_H
#define _PRINTLEVELS_H
const int print_none = 0;
const int print_low = 10;
const int print_medium = 20;
const int print_high = 30;
const int print_debug = 100;
const int print_default = print_medium;
const int print_never = 1000000;
#endif

56
src/util/stdio.h Normal file
View file

@ -0,0 +1,56 @@
#ifndef _USER_STDIO_H
#define _USER_STDIO_H
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// NAME
// stdio -- define logical units for standard I/O
//
// REVISION
// $Id$
//
// NOTES
// The common block must be initialized prior to using the I/O
// units. Currently the following points change these units:
//
// 1) Block data util_stdio_data [util_io.F] sets LuOut to 6 as a
// sensible default.
//
// 2) Function util_sgroup_set_ioname [util_sgroup.F] sets LuOut
// to a value based on the group number.
//
// 3) Function util_sgroup_unset_io [util_sgroup.F] closes LuOut.
//
// 4) Subroutine smd_group_set_io [smd_group.F] closes LuOut,
// resets it, and attaches it to a new file.
//
// 5) Subroutine smd_group_set_io_custom [smd_group.F] closes LuOut,
// resets it, and attaches it to a new file.
//
// 6) Subroutine smd_group_unset_io [smd_group.F] closes LuOut.
//
// This combination ensures that subgroup aware codes can arrange
// the I/O capabilities they need, while functionality that is
// not subgroup aware still works because of a proper default
// setting.
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
//
// This way we do not have to worry about the
// initialization/termination
//
char LuOut[511];
#endif
//
// A potentially useful tidbit: On Cray machines, units
// 100, 101, and 102 are always assigned to stdin, stdout, and
// stderr. They differ from 5, 6, and 0 in that they cannot be
// OPENed, and will not exist according to INQUIRE. Consequently,
// 100+ will _always_ correspond to the unix stdio streams regardless
// of what the application may do with 5/6/0
//
// Also note that on Crays, all of these units are _assigned_ but not
// preconnected. That means if you try to call something like flush
// on a unit that you have not written to previously (implicitly
// opeining it), it causes a fatal error.
#endif

11
src/util/util.h Normal file
View file

@ -0,0 +1,11 @@
#ifndef _UTIL_H
#define _UTIL_H
#include <stdbool.h>
#include "printlevels.h"
#include "util_maxlength.h"
const int nw_max_path_len = 255; // Maximum path len -> posix standard is what?
const int nw_max_path_len = MAXLENGTH; // Maximum path len -> posix standard is what?
#endif