pair_runner.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Feb 17, 05:58

pair_runner.cpp
View Options

	///////////////////////////////////////////////////////////////////////////////
	//
	// RuNNer pair style for LAMMPS
	// ----------------------------
	//
	// author: Andreas Singraber
	// date : 2015-03-10
	// email : andreas.singraber@univie.ac.at
	//
	// Copyright (c) 2013 - 2015 Andreas Singraber
	//
	////////////////////////////////////////////////////////////////////////////////

	#include <unistd.h>
	#include "math.h"
	#include "stdlib.h"
	#include "string.h"
	#include "pair_runner.h"
	#include "atom.h"
	#include "force.h"
	#include "domain.h"
	#include "comm.h"
	#include "neighbor.h"
	#include "neigh_list.h"
	#include "neigh_request.h"
	#include "memory.h"
	#include "error.h"
	#include "update.h"
	//#include "assert.h"

	using namespace LAMMPS_NS;

	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::PairRuNNer(LAMMPS *lmp) : Pair(lmp) {}

	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::~PairRuNNer()
	{
	if(allocated) {
	memory->destroy(setflag);
	memory->destroy(cutsq);
	}

	if(setup_completed == 1) {

	// L
	for(int i=0; i<ne; i++) {
	for(int j=0; j<nl+1; j++) {
	delete L[i][j];
	}
	}
	for(int i=0; i<ne; i++) {
	delete[] L[i];
	}
	delete[] L;

	// G
	for(int i=0; i<ne; i++) {
	for(int j=0; j<nsym[i]; j++) {
	delete G[i][j];
	}
	}
	for(int i=0; i<ne; i++) {
	delete[] G[i];
	}
	delete[] G;

	// GG
	for(int i=0; i<ne; i++) {
	for(int j=0; j<nsymg[i]; j++) {
	delete GG[i][j];
	}
	}
	for(int i=0; i<ne; i++) {
	delete[] GG[i];
	}
	delete[] GG;

	// a
	for(int i=0; i<na; i++) {
	delete a[i];
	}
	delete[] a;

	// e
	for(int i=0; i<ne; i++) {
	delete e[i];
	}
	delete[] e;

	// numnode, actfunc
	for(int i=0; i<ne; i++) delete[] numnode[i];
	delete[] numnode;
	for(int i=0; i<ne; i++) delete[] actfunc[i];
	delete[] actfunc;

	// nsym, nsymg
	delete[] nsym;
	delete[] nsymg;

	}

	}

	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::compute(int eflag, int vflag) {

	if(eflag \|\| vflag) ev_setup(eflag, vflag);
	else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;

	// allocate and set per-timestep arrays
	pt_allocset();

	// calculate neighbor list
	calc_neighbor_list();
	//calc_neighbor_list_lammps();

	// allocate space for force calculation
	init_forces();

	// calculate symmetry functions for all local atoms
	//calc_all_symfunc();
	calc_all_symfunc_group();

	// perform neural network calculation for all local atoms
	calc_nn();

	// calculate force contributions to all atoms
	calc_forces();

	// add free atom energies
	for(int i=0; i<nla; i++) {
	a[la[i]]->en += e[a[la[i]]->e]->atomic_energy;
	}

	// unit conversion for energies and forces
	for(int i=0; i<nla; i++) {
	a[la[i]]->en /= CONST_EV2HA;
	}
	for(int i=0; i<nla+atom->nghost; i++) {
	atom->f[i][0] *= CONST_A2B / CONST_EV2HA;
	atom->f[i][1] *= CONST_A2B / CONST_EV2HA;
	atom->f[i][2] *= CONST_A2B / CONST_EV2HA;
	}

	// sum up total energy
	entot = 0.0;
	for(int i=0; i<nla; i++) {
	entot += a[la[i]]->en;
	}

	// add total potential of this processor to global potential energy
	if(eflag_global) {
	ev_tally(0, 0, nla, 1, entot, 0.0, 0.0, 0.0, 0.0, 0.0);
	}

	if(eflag_atom) {
	for(int i=0; i<nla; i++) {
	eatom[i] = a[la[i]]->en;
	}
	}

	// if virial needed calculate via F dot r
	if(vflag_fdotr) {
	virial_fdotr_compute();
	}

	// print extrapolation warning summary if requested
	if(showewsum > 0 && update->ntimestep % showewsum == 0) {
	mpicollect_ew();
	if(comm->me == 0) {
	if(screen) {
	fprintf(screen, "# RuNNer EW summary: %ld %ld %f\n", update->ntimestep, ewcountsum, double(ewcountsum) / showewsum);
	}
	if(logfile) {
	fprintf(logfile, "# RuNNer EW summary: %ld %ld %f\n", update->ntimestep, ewcountsum, double(ewcountsum) / showewsum);
	}
	}
	ewcountsum = 0;
	}

	// delete per-timestep arrays
	pt_delete();

	}

	////////////////////////////////////////////////////////////////////////////////
	// allocate all arrays
	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::allocate() {

	int n = atom->ntypes;

	allocated = 1;
	memory->create(setflag, n+1, n+1, "pair:setflag");
	for(int i=1; i<=n; i++) {
	for(int j=i; j<=n; j++) {
	setflag[i][j] = 0;
	}
	}

	memory->create(cutsq, n+1, n+1, "pair:cutsq");

	setup_completed = 0;

	}

	////////////////////////////////////////////////////////////////////////////////
	// global settings
	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::settings(int narg, char **arg) {

	int iarg = 0;

	if(narg == 0) error->all(FLERR, "Illegal pair_style command");

	// default settings
	strcpy(dirRuNNer, "RuNNer");
	showew = 1;
	showewsum = 0;
	resetew = 0;
	maxew = 0;
	ntransformations = 0; //HELLSTROM


	while(iarg<narg) {
	// set RuNNer directory
	if(strcmp(arg[iarg], "dir") == 0) {
	if(iarg + 2 > narg) {
	error->all(FLERR, "Illegal pair_style command");
	}
	strcpy(dirRuNNer, arg[iarg+1]);
	iarg += 2;
	}
	// show extrapolation warnings
	else if(strcmp(arg[iarg], "showew") == 0) {
	if(iarg + 2 > narg) {
	error->all(FLERR, "Illegal pair_style command");
	}
	if(strcmp(arg[iarg+1], "yes") == 0) {
	showew = 1;
	}
	else if(strcmp(arg[iarg+1], "no") == 0) {
	showew = 0;
	}
	else {
	error->all(FLERR, "Illegal pair_style command");
	}
	iarg += 2;
	}
	// show extrapolation warning summary
	else if(strcmp(arg[iarg], "showewsum") == 0) {
	if(iarg + 2 > narg) {
	error->all(FLERR, "Illegal pair_style command");
	}
	sscanf(arg[iarg+1], "%d", &showewsum);
	iarg += 2;
	}
	// reset extrapolation warning counter
	else if(strcmp(arg[iarg], "resetew") == 0) {
	if(iarg + 2 > narg) {
	error->all(FLERR, "Illegal pair_style command");
	}
	if(strcmp(arg[iarg+1], "yes") == 0) {
	resetew = 1;
	}
	else if(strcmp(arg[iarg+1], "no") == 0) {
	resetew = 0;
	}
	else {
	error->all(FLERR, "Illegal pair_style command");
	}
	iarg += 2;
	}
	// maximum allowed extrapolation warnings
	else if(strcmp(arg[iarg], "maxew") == 0) {
	if(iarg + 2 > narg) {
	error->all(FLERR, "Illegal pair_style command");
	}
	sscanf(arg[iarg+1], "%d", &maxew);
	iarg += 2;
	}
	else if(strcmp(arg[iarg], "transform") == 0) { //HELLSTROM
	// transform 4 1 ===> atoms of type 4 will be "chemically" treated as atoms of type 1 (symmetry functions etc. are calculated as for atom type 1)
	// the benefit is that you can put more atom types in the structure file if you want several isotopes of the same element
	// it is necessary to use a different atom type because LAMMPS only allows masses per atom type (not per individual atom, other than for some specific strange atom styles)
	// the user does not need to change anything in input.nn, weights.xxx.data, scaling.data
	if(iarg + 3 > narg) {
	error->all(FLERR, "Illegal pair_style (transform) command");
	}
	ntransformations += 1;
	if (ntransformations > MAXTRANSFORMATIONS) error->all(FLERR, "Too many transformations!");
	sscanf(arg[iarg+1], "%d", &transformfrom[ntransformations-1]);
	sscanf(arg[iarg+2], "%d", &transformto[ntransformations-1]);
	iarg += 3;
	}
	else {
	error->all(FLERR, "Illegal pair_style command");
	}
	}
	sprintf(dirRuNNer, "%s/", dirRuNNer);

	// show settings
	if(comm->me == 0) {
	if(screen) {
	fprintf(screen, "RuNNer pair style directory : %s\n", dirRuNNer);
	if (showew == 1) fprintf(screen, "Show extrapolation warnings : %s\n", "yes");
	else if(showew == 0) fprintf(screen, "Show extrapolation warnings : %s\n", "no");
	if (showewsum > 0) fprintf(screen, "Show extrapolation warning summary : every %d timesteps\n", showewsum);
	else fprintf(screen, "Show extrapolation warning summary : %s\n", "no");
	if (resetew == 1) fprintf(screen, "Reset extrapolation warning counter : %s\n", "yes");
	else if(resetew == 0) fprintf(screen, "Reset extrapolation warning counter : %s\n", "no");
	fprintf(screen, "Maximum allowed extrapolation warnings: %d\n", maxew);
	}
	if(logfile) {
	fprintf(logfile, "RuNNer pair style directory : %s\n", dirRuNNer);
	if (showew == 1) fprintf(logfile, "Show extrapolation warnings : %s\n", "yes");
	else if(showew == 0) fprintf(logfile, "Show extrapolation warnings : %s\n", "no");
	if (showewsum > 0) fprintf(logfile, "Show extrapolation warning summary : every %d timesteps\n", showewsum);
	else fprintf(logfile, "Show extrapolation warning summary : %s\n", "no");
	if (resetew == 1) fprintf(logfile, "Reset extrapolation warning counter : %s\n", "yes");
	else if(resetew == 0) fprintf(logfile, "Reset extrapolation warning counter : %s\n", "no");
	fprintf(logfile, "Maximum allowed extrapolation warnings: %d\n", maxew);
	}
	}

	}

	////////////////////////////////////////////////////////////////////////////////
	// set coeffs for one or more type pairs
	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::coeff(int narg, char **arg) {

	if(narg > 3) error->all(FLERR, "Incorrect args for pair coefficients");
	if(!allocated) allocate();

	int ilo,ihi,jlo,jhi;
	force->bounds(FLERR, arg[0], atom->ntypes, ilo, ihi);
	force->bounds(FLERR, arg[1], atom->ntypes, jlo, jhi);

	cut_global = force->numeric(FLERR, arg[2]);

	int count = 0;
	for(int i=ilo; i<=ihi; i++) {
	for(int j=MAX(jlo,i); j<=jhi; j++) {
	setflag[i][j] = 1;
	count++;
	}
	}

	if(count == 0) error->all(FLERR, "Incorrect args for pair coefficients");

	}

	////////////////////////////////////////////////////////////////////////////////
	// initialization of pair style
	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::init_style() {

	int irequest = neighbor->request((void *) this);
	neighbor->requests[irequest]->pair = 1;
	neighbor->requests[irequest]->half = 0;
	neighbor->requests[irequest]->full = 1;

	if(setup_completed == 0) {
	if(comm->me == 0) {
	if(screen) {
	fprintf(screen, "-----------------------------------------------------------------------\n");
	fprintf(screen, "Starting RuNNer pair style setup...\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	fprintf(logfile, "Starting RuNNer pair style setup...\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}

	// INITIALIZE SOME VARIABLES
	init_vars();

	// GET INFORMATION FROM input.nn
	if(comm->me == 0) {
	analyze_input_nn();
	}

	// MPI: DISTRIBUTE NN SETUP
	mpidistribute_nn();

	// ALLOCATE SPACE FOR NETWORK LAYERS
	init_layers();

	// READ NEURAL NETWORK WEIGHTS FROM weights.???.data
	if(comm->me == 0) {
	read_weights();
	}

	// MPI: DISTRIBUTE NEURAL NETWORK WEIGHTS
	mpidistribute_nnweights();

	// SORT SYMMETRY FUNCTIONS
	if(comm->me == 0) {
	sort_symfuncs();
	}

	// READ IN SCALING VALUES FROM scaling.data
	if(comm->me == 0) {
	read_scaling();
	}

	// MPI: DISTRIBUTE SYMMETRY FUNCTIONS
	mpidistribute_symfunc();

	// FIND MAXIMUM GLOBAL CUTOFF
	rcmax = rc_max();

	// CREATE SF GROUPS
	create_symfunc_groups();

	// ALLOCATE RuNNer_atom ARRAY
	na = atom->natoms;
	a = new RuNNer_atom*[na];
	for(int i=0; i<na; i++) a[i] = new RuNNer_atom(*this);

	setup_completed = 1;

	if(comm->me == 0) {
	if(screen) {
	fprintf(screen, "RuNNer pair style setup completed.\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "RuNNer pair style setup completed.\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	}

	}

	////////////////////////////////////////////////////////////////////////////////
	// init for one type pair i,j and corresponding j,i
	////////////////////////////////////////////////////////////////////////////////

	double PairRuNNer::init_one(int i, int j) {

	if(setflag[i][j] == 0) {
	}

	if(offset_flag) {
	}
	else {
	}

	return cut_global;

	}

	////////////////////////////////////////////////////////////////////////////////
	// proc 0 writes to restart file
	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::write_restart(FILE *fp) {

	write_restart_settings(fp);

	int i,j;
	for(i=1; i<=atom->ntypes; i++)
	for(j=i; j<=atom->ntypes; j++) {
	fwrite(&setflag[i][j], sizeof(int), 1, fp);
	if(setflag[i][j]) {
	}
	}

	}

	////////////////////////////////////////////////////////////////////////////////
	// proc 0 reads from restart file, bcasts
	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::read_restart(FILE *fp) {

	read_restart_settings(fp);

	allocate();

	int i, j;
	int me = comm->me;
	for(i=1; i<=atom->ntypes; i++) {
	for(j=i; j<=atom->ntypes; j++) {
	if(me == 0) fread(&setflag[i][j], sizeof(int), 1, fp);
	MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
	if(setflag[i][j]) {
	if(me == 0) {
	}
	}
	}
	}

	}

	////////////////////////////////////////////////////////////////////////////////
	// proc 0 writes to restart file
	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::write_restart_settings(FILE *fp) {

	fwrite(&offset_flag, sizeof(int), 1, fp);
	fwrite(&mix_flag, sizeof(int), 1, fp);

	}

	////////////////////////////////////////////////////////////////////////////////
	// proc 0 reads from restart file, bcasts
	////////////////////////////////////////////////////////////////////////////////

	void PairRuNNer::read_restart_settings(FILE *fp) {

	if(comm->me == 0) {
	fread(&offset_flag, sizeof(int), 1, fp);
	fread(&mix_flag, sizeof(int), 1, fp);
	}
	MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
	MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);

	}

	////////////////////////////////////////////////////////////////////////////////

	double PairRuNNer::single(int i, int j, int itype, int jtype, double rsq,
	double factor_coul, double factor_lj, double &fforce) {

	error->all(FLERR, "PairRuNNer::single not possible for pair_style runner");

	return 0.0;

	}


	////////////////////////////////////////////////////////////////////////////////
	// PairRuNNer member functions
	////////////////////////////////////////////////////////////////////////////////

	// initialize some variables
	void PairRuNNer::init_vars() {

	setup_completed = 0;
	na = 0;
	nla = 0;
	ne = 0;
	nl = 0;
	scalesym = 0;
	centersym = 0;
	normnodes = 0;
	ewcount = 0;
	ewcountsum = 0;
	rcmax = 0.0;
	entot = 0.0;
	scalemin = 0.0;
	scalemax = 0.0;

	}

	int PairRuNNer::transform_atom_type(int original_atom_type) const { //HELLSTROM
	for (int j=0; j<ntransformations; j++) {
	if (transformfrom[j] == original_atom_type) {
	//fprintf(stderr, "Transforming atom from %d to %d\n", original_atom_type, transformto[j]);
	return transformto[j];
	}
	}
	return original_atom_type;
	}

	// allocate and set per-timestep arrays
	void PairRuNNer::pt_allocset() {

	int acount = 0;
	int index = 0;

	// create list of local atoms
	nla = atom->nlocal;
	la = new int[nla];
	nlape = new int[ne];
	for(int i=0; i<ne; i++) nlape[i] = 0;
	for(int i=0; i<nla; i++) {
	la[i] = atom->tag[i] - 1;
	//nlape[atom->type[i]-1]++;
	nlape[transform_atom_type(atom->type[i])-1]++; //HELLSTROM
	}

	// allocate RuNNer_neighbor array
	N = new RuNNer_neighbor*[nla];
	for(int i=0; i<nla; i++) N[i] = new RuNNer_neighbor(*this);

	// transfer local atom posisitons to RuNNer_atom array
	for(int i=0; i<nla; i++) {
	index = la[i];
	a[index]->x = CONST_A2B * atom->x[i][0];
	a[index]->y = CONST_A2B * atom->x[i][1];
	a[index]->z = CONST_A2B * atom->x[i][2];
	//a[index]->e = atom->type[i]-1;
	a[index]->e = transform_atom_type(atom->type[i])-1; //HELLSTROM
	}

	// allocate SF value array
	laipe = new int*[ne];
	G_svalue = new double**[ne];
	for(int i=0; i<ne; i++) {
	laipe[i] = new int[nlape[i]];
	G_svalue[i] = new double*[nlape[i]];
	for(int j=0; j<nlape[i]; j++) {
	laipe[i][j] = 0;
	G_svalue[i][j] = new double[nsym[i]];
	for(int k=0; k<nsym[i]; k++) {
	G_svalue[i][j][k] = 0.0;
	}
	}
	}

	// set per-element atom index array
	for(int i=0; i<ne; i++) {
	acount = 0;
	for(int j=0; j<nla; j++) {
	if(a[la[j]]->e == i) {
	laipe[i][acount] = j;
	acount++;
	}
	}
	}

	// reset extrapolation warning counter if requested
	if(resetew == 1) {
	ewcount = 0;
	}

	}

	// delete per-timestep arrays
	void PairRuNNer::pt_delete() {

	// free space
	for(int i=0; i<ne; i++) {
	for(int j=0; j<nlape[i]; j++) {
	delete[] G_svalue[i][j];
	}
	delete[] G_svalue[i];
	delete[] laipe[i];
	}
	delete[] G_svalue;
	delete[] laipe;
	delete[] la;
	delete[] nlape;
	for(int i=0; i<nla; i++) {
	delete N[i];
	delete F[i];
	}
	delete[] N;
	delete[] F;

	}

	// sort elements according to nuclear charge
	void PairRuNNer::sort_elements() {

	int i;
	RuNNer_element temp(*this);

	i = 0;
	do {
	if(e[i]->nucch > e[i+1]->nucch) {
	temp = *e[i];
	e[i] = e[i+1];
	*e[i+1] = temp;
	i = -1;
	}
	i++;
	} while(i<ne-1);

	return;

	}

	// find element number for a given element string
	int PairRuNNer::find_element(char *tmpname) {

	int i;

	for(i=0; i<ne; i++) {
	if(strcmp(e[i]->name, tmpname) == 0) return i;
	}
	fprintf(stderr, "ERROR: element %s not specified!\n", tmpname);
	exit(EXIT_FAILURE);

	}

	// initialize force classes
	void PairRuNNer::init_forces() {

	F = new RuNNer_forces*[nla];
	for(int i=0; i<nla; i++) {
	F[i] = new RuNNer_forces(*this, nl, numnode[N[i]->eatom], N[i]->eatom, i, N[i]->num);
	}

	return;

	}

	// calculate forces on local and ghost atoms from saved data
	void PairRuNNer::calc_forces() {

	int i = 0;
	int j = 0;
	int ca = 0;
	int ica = 0;
	int neigh = 0;
	int eneigh = 0;
	int ineighcaind = 0;

	// loop over all local + ghost atoms
	for(ica=0; ica<nla+atom->nghost; ica++) {
	// loop over all possible neighbors = local atoms
	for(i=0; i<nla; i++) {
	neigh = la[i];
	eneigh = a[neigh]->e;
	ineighcaind = N[i]->find_nindex(ica);
	if(ineighcaind >= 0) {
	// if(N[i]->ilal2.find(ica)!=N[i]->ilal2.end()) {
	// ineighcaind = N[i]->ilal2.at(ica);
	// //assert(ineighcaind==N[i]->find_nindex(ica));
	for(j=0; j<nsym[eneigh]; j++) {
	atom->f[ica][0] -= F[i]->dEdG[j] * F[i]->dGdx[j][ineighcaind] * G[eneigh][j]->sf;
	atom->f[ica][1] -= F[i]->dEdG[j] * F[i]->dGdy[j][ineighcaind] * G[eneigh][j]->sf;
	atom->f[ica][2] -= F[i]->dEdG[j] * F[i]->dGdz[j][ineighcaind] * G[eneigh][j]->sf;
	}
	}
	}
	}
	// add force contribution of central atom itself
	for(ica=0; ica<nla; ica++) {
	ca = la[ica];
	for(j=0; j<nsym[a[ca]->e]; j++) {
	atom->f[ica][0] -= F[ica]->dEdG[j] * F[ica]->dGdx[j][N[ica]->num] * G[a[ca]->e][j]->sf;
	atom->f[ica][1] -= F[ica]->dEdG[j] * F[ica]->dGdy[j][N[ica]->num] * G[a[ca]->e][j]->sf;
	atom->f[ica][2] -= F[ica]->dEdG[j] * F[ica]->dGdz[j][N[ica]->num] * G[a[ca]->e][j]->sf;
	}
	}

	return;

	}

	// allocate space for layers
	void PairRuNNer::init_layers() {

	L = new RuNNer_layer**[ne];
	for(int i=0; i<ne; i++) {
	L[i] = new RuNNer_layer*[nl+1];
	}
	for(int i=0; i<ne; i++) {
	for(int j=0; j<nl+1; j++) {
	L[i][j] = new RuNNer_layer(*this, numnode[i][j], numnode[i][j+1], actfunc[i][j+1]);
	if(comm->me==0) {
	if(screen) {
	fprintf(screen, "Element %2d (%2s): layer %2d allocated, dim: (%3d,%3d), act. func.: %2d\n",
	i, e[i]->name, j, numnode[i][j], numnode[i][j+1], actfunc[i][j+1]);
	}
	if(logfile) {
	fprintf(logfile, "Element %2d (%2s): layer %2d allocated, dim: (%3d,%3d), act. func.: %2d\n",
	i, e[i]->name, j, numnode[i][j], numnode[i][j+1], actfunc[i][j+1]);
	}
	}
	}
	}
	if(comm->me==0) {
	if(screen) fprintf(screen, "-----------------------------------------------------------------------\n");
	if(logfile) fprintf(logfile, "-----------------------------------------------------------------------\n");
	}

	return;

	}

	// open weights.???.data file and call get_weights
	void PairRuNNer::read_weights() {

	FILE *fp;
	char fname[LLINE+16];
	int nweight = 0;
	int nbias = 0;

	for(int i=0; i<ne; i++) {
	sprintf(fname, "%sweights.%03d.data", dirRuNNer, e[i]->nucch);
	if(screen) fprintf(screen, "Element %2d (%2s): Reading %s...\n", i, e[i]->name, fname);
	if(logfile) fprintf(logfile, "Element %2d (%2s): Reading %s...\n", i, e[i]->name, fname);
	fp = fopen(fname, "r");
	if( fp == NULL ) {
	fprintf(stderr, "ERROR: Could not open %s.\n", fname);
	exit(EXIT_FAILURE);
	}
	else {
	nweight = 0;
	nbias = 0;
	for(int j=0; j<nl+1; j++) {
	L[i][j]->get_weights(fp);
	nweight += numnode[i][j] * numnode[i][j+1];
	nbias += numnode[i][j+1];
	}
	fclose(fp);
	}
	if(screen) {
	fprintf(screen, "Element %2d (%2s): weights: %6d, bias: %4d, total: %6d\n",
	i, e[i]->name, nweight, nbias, nweight+nbias);
	}
	if(logfile) {
	fprintf(logfile, "Element %2d (%2s): weights: %6d, bias: %4d, total: %6d\n",
	i, e[i]->name, nweight, nbias, nweight+nbias);
	}
	}
	if(screen) fprintf(screen, "-----------------------------------------------------------------------\n");
	if(logfile) fprintf(logfile, "-----------------------------------------------------------------------\n");

	return;

	}

	// calculate NN for all atoms
	void PairRuNNer::calc_nn() {

	int ca = 0;

	// loop over all elements
	for(int i=0; i<ne; i++) {
	// loop over all atoms of element i
	for(int ica=0; ica<nlape[i]; ica++) {
	ca = laipe[i][ica];
	// set symmetry functions values as input layer nodes
	for(int isym=0; isym<nsym[i]; isym++) {
	L[i][0]->input[isym] = G_svalue[i][ica][isym];
	}
	// calculate all nodes of NN
	L[i][0]->calc_all_nodes();
	for(int k=1;k<nl+1;k++) {
	L[i][k]->get_input(L[i][k-1]);
	L[i][k]->calc_all_nodes();
	}
	// save data needed for force calculation
	F[ca]->get_dfdx(L[i]);
	F[ca]->calc_dEdG(L[i]);
	// atomic energy contribution is value of output node
	a[la[ca]]->en = L[i][nl]->node[0];
	}
	}

	return;

	}

	// calculate neighbor list for all atoms
	void PairRuNNer::calc_neighbor_list() {

	int i = 0;

	// loop over all atoms
	//#ifdef _OPENMP
	// #pragma omp parallel for private(i)
	//#endif
	for(i=0; i<nla; i++) {
	N[i]->calc(a, la[i]);
	}

	return;

	}

	// calculate neighbor list for all atoms
	void PairRuNNer::calc_neighbor_list_lammps() {

	int i = 0;
	int inum = list->inum;
	int* ilist = list->ilist;

	// loop over all atoms
	for(int ii=0; ii<inum; ii++) {
	i = ilist[ii];
	N[i]->calcl(a, ii);
	}

	return;

	}

	// Read necessary information from input.nn
	void PairRuNNer::analyze_input_nn() {

	FILE *fp;
	int i, j;
	int *isf;
	int tmpe;
	int tmptype;
	char fname[LLINE+9];
	char line[LLINE];
	char gsfline[LLINE];
	char keyword[LLINE];
	char tmpename[3];
	int ierr;
	double tmpd;

	ne = 0;
	nl = 0;
	cutoff_type = 1;
	sprintf(fname, "%sinput.nn", dirRuNNer);
	fp = fopen(fname, "r");
	if(screen) {
	fprintf(screen, "Analyzing input.nn...\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "Analyzing input.nn...\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	if(fp == NULL) {
	fprintf(stderr, "ERROR: Could not open input.nn.\n");
	exit(EXIT_FAILURE);
	}
	else {
	while(fgets(line, LLINE, fp) != NULL) {
	keyword[0]='\0';
	ierr = sscanf(line, "%s", keyword);
	// detect number of elements
	if( (ne == 0) && (strcmp(keyword, "number_of_elements") == 0) ) {
	ierr = sscanf(line, "%*s%d", &ne);
	e = new RuNNer_element*[ne];
	nsym = new int[ne];
	for(i=0; i<ne; i++) {
	e[i] = new RuNNer_element(*this);
	nsym[i] = 0;
	}
	}
	// read names of elements
	if( (ne != 0) && (strcmp(keyword, "elements") == 0) ) {
	ierr = sscanf(line, "%*s %[^\n]", line);
	for(i=0; i<ne; i++) {
	ierr = sscanf(line, "%s %[^\n]", e[i]->name, line);
	e[i]->getnucch();
	}
	if( ne > 1) sort_elements();
	if(screen) {
	fprintf(screen, "Number of elements: %2d\n", ne);
	for(i=0; i<ne; i++) fprintf(screen, "Element %2d: %2s (%2d)\n", i, e[i]->name, e[i]->nucch);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "Number of elements: %2d\n", ne);
	for(i=0; i<ne; i++) fprintf(logfile, "Element %2d: %2s (%2d)\n", i, e[i]->name, e[i]->nucch);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	// read single symmetry functions
	if( (ne != 0) && (strcmp(keyword, "symfunction_short") == 0) ) {
	strcpy(tmpename, "?");
	ierr = sscanf(line, "%*s%s", tmpename);
	nsym[find_element(tmpename)]++;
	}
	// read element symmetry functions
	if( (ne != 0) && (strcmp(keyword, "element_symfunction_short") == 0) ) {
	tmpe = 0;
	strcpy(tmpename, "?");
	ierr = sscanf(line, "%*s%s%d", tmpename, &tmptype);
	tmpe = find_element(tmpename);
	if(tmptype == 2) {
	nsym[tmpe] += ne;
	}
	else if(tmptype == 3) {
	nsym[tmpe] += ne * (ne + 1) / 2;
	}
	else {
	fprintf(stderr, "ERROR: symmetry function type %d not implemented.\n", tmptype);
	exit(EXIT_FAILURE);
	}
	}
	// read global symmetry functions
	if( (ne != 0) && (strcmp(keyword, "global_symfunction_short") == 0) ) {
	ierr = sscanf(line, "%*s%d", &tmptype);
	if(tmptype == 2) {
	for(int i=0; i<ne; i++) {
	nsym[i] += ne;
	}
	}
	else if(tmptype == 3) {
	for(int i=0; i<ne; i++) {
	nsym[i] += ne * (ne + 1) / 2;
	}
	}
	else {
	fprintf(stderr, "ERROR: symmetry function type %d not implemented.\n", tmptype);
	exit(EXIT_FAILURE);
	}
	}
	// detect number of layers
	if( (ne != 0) && (strcmp(keyword, "global_hidden_layers_short") == 0) ) {
	ierr = sscanf(line, "%*s%d", &nl);
	if(screen) {
	fprintf(screen, "Number of hidden layers: %2d\n", nl);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "Number of hidden layers: %2d\n", nl);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	numnode = new int*[ne];
	actfunc = new int*[ne];
	for(i=0; i<ne; i++) {
	numnode[i] = new int[nl+2];
	actfunc[i] = new int[nl+2];
	for(j=0; j<nl+2; j++) {
	numnode[i][j]=0;
	actfunc[i][j]=0;
	}
	numnode[i][nl+1]=1;
	actfunc[i][0]=0;
	}
	}
	// detect activation function for each layer
	if( (nl != 0) && (strcmp(keyword, "global_nodes_short") == 0) ) {
	ierr = sscanf(line, "%*s %[^\n]", line);
	for(i=1; i<nl; i++) {
	ierr = sscanf(line, "%d %[^\n]", &numnode[0][i], line);
	for(j=1; j<ne; j++) numnode[j][i] = numnode[0][i];
	}
	ierr = sscanf(line, "%d", &numnode[0][nl]);
	for(j=1; j<ne; j++) numnode[j][nl] = numnode[0][nl];
	for(i=0; i<ne; i++) {
	for(j=1; j<nl+2; j++) {
	if(screen) fprintf(screen, "Element %2d (%2s): layer %2d: %3d nodes\n", i, e[i]->name, j-1, numnode[i][j]);
	if(logfile) fprintf(logfile, "Element %2d (%2s): layer %2d: %3d nodes\n", i, e[i]->name, j-1, numnode[i][j]);
	}
	}
	}
	// detect number of nodes of each layer
	if( (nl != 0) && (strcmp(keyword, "global_activation_short") == 0) ) {
	char tmpact[2];
	ierr = sscanf(line, "%*s %[^\n]", line);
	for(i=1; i<nl+1; i++) {
	ierr = sscanf(line, "%s %[^\n]", tmpact, line);
	if(strcmp(tmpact, "l") == 0) actfunc[0][i]=1;
	if(strcmp(tmpact, "t") == 0) actfunc[0][i]=2;
	for(j=1; j<ne; j++) actfunc[j][i] = actfunc[0][i];
	}
	ierr = sscanf(line, "%s", tmpact);
	if(strcmp(tmpact, "l") == 0) actfunc[0][nl+1] = 1;
	if(strcmp(tmpact, "t") == 0) actfunc[0][nl+1] = 2;
	for(j=1; j<ne; j++) actfunc[j][nl+1] = actfunc[0][nl+1];
	for(i=0; i<ne; i++) {
	for(j=1; j<nl+2; j++) {
	if(screen) {
	fprintf(screen, "Element %2d (%2s): layer %2d: activation function: %2d\n",
	i, e[i]->name, j-1, actfunc[i][j]);
	}
	if(logfile) {
	fprintf(logfile, "Element %2d (%2s): layer %2d: activation function: %2d\n",
	i, e[i]->name, j-1, actfunc[i][j]);
	}
	}
	}
	if(screen) fprintf(screen, "-----------------------------------------------------------------------\n");
	if(logfile) fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	// detect scale_symmetry_functions
	if( strcmp(keyword, "scale_symmetry_functions") == 0 ) {
	scalesym = 1;
	if(screen) {
	fprintf(screen, "scale_symmetry_functions detected.\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "scale_symmetry_functions detected.\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	// detect center_symmetry_functions
	if( strcmp(keyword, "center_symmetry_functions") == 0 ) {
	centersym = 1;
	if(screen) {
	fprintf(screen, "center_symmetry_functions detected.\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "center_symmetry_functions detected.\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	// set scalemin value
	if( strcmp(keyword, "scale_min_short") == 0 ) {
	ierr = sscanf(line, "%*s%lf", &scalemin);
	if(screen) {
	fprintf(screen, "scalemin = %f\n", scalemin);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "scalemin = %f\n", scalemin);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	// set scalemax value
	if( strcmp(keyword, "scale_max_short") == 0 ) {
	ierr = sscanf(line, "%*s%lf", &scalemax);
	if(screen) {
	fprintf(screen, "scalemax = %f\n", scalemax);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "scalemax = %f\n", scalemax);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	// detect normalize_nodes
	if( strcmp(keyword, "normalize_nodes") == 0 ) {
	normnodes = 1;
	if(screen) {
	fprintf(screen, "normalize_nodes detected.\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "normalize_nodes detected.\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	// detect cutoff type
	if( (strcmp(keyword, "cutoff_type") == 0) ) {
	ierr = sscanf(line, "%*s%d", &cutoff_type);
	if(screen) {
	fprintf(screen, "cutoff type = %d\n", cutoff_type);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "cutoff type = %d\n", cutoff_type);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	// read atomic energies
	if( (strcmp(keyword, "atom_energy") == 0) ) {
	ierr = sscanf(line, "%*s%s%lf", tmpename, &tmpd);
	tmpe = find_element(tmpename);
	e[tmpe]->atomic_energy = tmpd;
	if(screen) {
	fprintf(screen, "Free atom energy of element %2d (%2s) = %16.8E\n", tmpe, e[tmpe]->name, e[tmpe]->atomic_energy);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "Free atom energy of element %2d (%2s) = %16.8E\n", tmpe, e[tmpe]->name, e[tmpe]->atomic_energy);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	}
	if(screen) {
	for(i=0; i<ne; i++) fprintf(screen, "Element %2d (%2s): %d symmetry functions\n", i, e[i]->name, nsym[i]);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	for(i=0; i<ne; i++) fprintf(logfile, "Element %2d (%2s): %d symmetry functions\n", i, e[i]->name, nsym[i]);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	G = new RuNNer_symfunc**[ne];
	for(i=0; i<ne; i++) {
	G[i] = new RuNNer_symfunc*[nsym[i]];
	for(j=0; j<nsym[i]; j++) G[i][j] = new RuNNer_symfunc(*this);
	numnode[i][0] = nsym[i];
	}
	// read input.nn again and set symmetry function parameters
	rewind(fp);
	isf = new int[ne];
	for(i=0; i<ne; i++) isf[i] = 0;
	while(fgets(line, LLINE, fp) != NULL) {
	keyword[0] = '\0';
	ierr = sscanf(line, "%s", keyword);
	if( (ne != 0) && (strcmp(keyword, "symfunction_short") == 0) ) {
	tmpe = 0;
	strcpy(tmpename, "?");
	ierr = sscanf(line, "%*s%s", tmpename);
	tmpe = find_element(tmpename);
	ierr = sscanf(line, "%*s %[^\n]", line);
	G[tmpe][isf[tmpe]]->setup(line);
	isf[tmpe]++;
	}
	if( (ne != 0) && (strcmp(keyword, "element_symfunction_short") == 0) ) {
	tmpe = 0;
	strcpy(tmpename, "?");
	ierr = sscanf(line, "%*s%s%d", tmpename, &tmptype);
	tmpe = find_element(tmpename);
	ierr = sscanf(line, "%s %s %*d %[^\n]", line);
	if(tmptype == 2) {
	for(int i=0; i<ne; i++) {
	sprintf(gsfline, "%s %d %s %s", tmpename, tmptype, e[i]->name, line);
	G[tmpe][isf[tmpe]]->setup(gsfline);
	isf[tmpe]++;
	}
	}
	else if(tmptype == 3) {
	for(int i=0; i<ne; i++) {
	for(int j=i; j<ne; j++) {
	sprintf(gsfline, "%s %d %s %s %s", tmpename, tmptype, e[i]->name, e[j]->name, line);
	G[tmpe][isf[tmpe]]->setup(gsfline);
	isf[tmpe]++;
	}
	}
	}
	}
	if( (ne != 0) && (strcmp(keyword, "global_symfunction_short") == 0) ) {
	tmpe = 0;
	ierr = sscanf(line, "%*s%d", &tmptype);
	ierr = sscanf(line, "%s %d %[^\n]", line);
	if(tmptype == 2) {
	for(int i=0; i<ne; i++) {
	for(int j=0; j<ne; j++) {
	sprintf(gsfline, "%s %d %s %s", e[i]->name, tmptype, e[j]->name, line);
	G[i][isf[i]]->setup(gsfline);
	isf[i]++;
	}
	}
	}
	else if(tmptype == 3) {
	for(int i=0; i<ne; i++) {
	for(int j=0; j<ne; j++) {
	for(int k=j; k<ne; k++) {
	sprintf(gsfline, "%s %d %s %s %s", e[i]->name, tmptype, e[j]->name, e[k]->name, line);
	G[i][isf[i]]->setup(gsfline);
	isf[i]++;
	}
	}
	}
	}
	}
	}
	delete[] isf;
	fclose(fp);
	}
	if(screen) {
	fprintf(screen, "Finished analyzing input.nn.\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "Finished analyzing input.nn.\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}

	return;

	}

	// sort symmetry functions
	void PairRuNNer::sort_symfuncs() {

	int i, j;
	RuNNer_symfunc temp(*this);

	if(screen) fprintf(screen, "Sorting according to type.\n");
	if(logfile) fprintf(logfile, "Sorting according to type.\n");
	for(i=0; i<ne; i++) {
	if(nsym[i]>1) {
	j = 0;
	do {
	if(G[i][j]->type > G[i][j+1]->type) {
	temp = *G[i][j];
	G[i][j] = G[i][j+1];
	*G[i][j+1] = temp;
	j = -1;
	}
	j++;
	} while(j<nsym[i]-1);
	}
	}

	if(screen) fprintf(screen, "Sorting according to rc.\n");
	if(logfile) fprintf(logfile, "Sorting according to rc.\n");
	for(i=0; i<ne; i++) {
	if(nsym[i]>1) {
	j = 0;
	do {
	if( (G[i][j]->type == G[i][j+1]->type) &&
	(G[i][j]->rc > G[i][j+1]->rc) ) {
	temp = *G[i][j];
	G[i][j] = G[i][j+1];
	*G[i][j+1] = temp;
	j = -1;
	}
	j++;
	} while(j<nsym[i]-1);
	}
	}

	if(screen) fprintf(screen, "Sorting according to eta.\n");
	if(logfile) fprintf(logfile, "Sorting according to eta.\n");
	for(i=0; i<ne; i++) {
	if(nsym[i]>1) {
	j = 0;
	do {
	if( (G[i][j]->type == G[i][j+1]->type) &&
	(G[i][j]->rc == G[i][j+1]->rc) &&
	(G[i][j]->eta > G[i][j+1]->eta) ) {
	temp = *G[i][j];
	G[i][j] = G[i][j+1];
	*G[i][j+1] = temp;
	j = -1;
	}
	j++;
	} while(j<nsym[i]-1);
	}
	}

	if(screen) fprintf(screen, "Sorting according to zeta.\n");
	if(logfile) fprintf(logfile, "Sorting according to zeta.\n");
	for(i=0; i<ne; i++) {
	if(nsym[i]>1) {
	j = 0;
	do {
	if( (G[i][j]->type == G[i][j+1]->type) &&
	(G[i][j]->rc == G[i][j+1]->rc) &&
	(G[i][j]->eta == G[i][j+1]->eta) &&
	(G[i][j]->zeta > G[i][j+1]->zeta) ) {
	temp = *G[i][j];
	G[i][j] = G[i][j+1];
	*G[i][j+1] = temp;
	j = -1;
	}
	j++;
	} while(j<nsym[i]-1);
	}
	}

	if(screen) fprintf(screen, "Sorting according to lambda.\n");
	if(logfile) fprintf(logfile, "Sorting according to lambda.\n");
	for(i=0; i<ne; i++) {
	if(nsym[i]>1) {
	j = 0;
	do {
	if( (G[i][j]->type == G[i][j+1]->type) &&
	(G[i][j]->rc == G[i][j+1]->rc) &&
	(G[i][j]->eta == G[i][j+1]->eta) &&
	(G[i][j]->zeta == G[i][j+1]->zeta) &&
	(G[i][j]->lambda > G[i][j+1]->lambda) ) {
	temp = *G[i][j];
	G[i][j] = G[i][j+1];
	*G[i][j+1] = temp;
	j = -1;
	}
	j++;
	} while(j<nsym[i]-1);
	}
	}

	if(screen) fprintf(screen, "Sorting according to rs.\n");
	if(logfile) fprintf(logfile, "Sorting according to rs.\n");
	for(i=0; i<ne; i++) {
	if(nsym[i]>1) {
	j = 0;
	do {
	if( (G[i][j]->type == G[i][j+1]->type) &&
	(G[i][j]->rc == G[i][j+1]->rc) &&
	(G[i][j]->eta == G[i][j+1]->eta) &&
	(G[i][j]->zeta == G[i][j+1]->zeta) &&
	(G[i][j]->lambda == G[i][j+1]->lambda) &&
	(G[i][j]->rs > G[i][j+1]->rs) ) {
	temp = *G[i][j];
	G[i][j] = G[i][j+1];
	*G[i][j+1] = temp;
	j = -1;
	}
	j++;
	} while(j<nsym[i]-1);
	}
	}

	if(screen) fprintf(screen, "Sorting according to e1.\n");
	if(logfile) fprintf(logfile, "Sorting according to e1.\n");
	for(i=0; i<ne; i++) {
	if(nsym[i]>1) {
	j = 0;
	do {
	if( (G[i][j]->type == G[i][j+1]->type) &&
	(G[i][j]->rc == G[i][j+1]->rc) &&
	(G[i][j]->eta == G[i][j+1]->eta) &&
	(G[i][j]->zeta == G[i][j+1]->zeta) &&
	(G[i][j]->lambda == G[i][j+1]->lambda) &&
	(G[i][j]->rs == G[i][j+1]->rs) &&
	(G[i][j]->e1 > G[i][j+1]->e1) ) {
	temp = *G[i][j];
	G[i][j] = G[i][j+1];
	*G[i][j+1] = temp;
	j = -1;
	}
	j++;
	} while(j<nsym[i]-1);
	}
	}

	if(screen) fprintf(screen, "Sorting according to e2.\n");
	if(logfile) fprintf(logfile, "Sorting according to e2.\n");
	for(i=0; i<ne; i++) {
	if(nsym[i]>1) {
	j = 0;
	do {
	if( (G[i][j]->type == G[i][j+1]->type) &&
	(G[i][j]->rc == G[i][j+1]->rc) &&
	(G[i][j]->eta == G[i][j+1]->eta) &&
	(G[i][j]->zeta == G[i][j+1]->zeta) &&
	(G[i][j]->lambda == G[i][j+1]->lambda) &&
	(G[i][j]->rs == G[i][j+1]->rs) &&
	(G[i][j]->e1 == G[i][j+1]->e1) &&
	(G[i][j]->e2 > G[i][j+1]->e2) ) {
	temp = *G[i][j];
	G[i][j] = G[i][j+1];
	*G[i][j+1] = temp;
	j = -1;
	}
	j++;
	} while(j<nsym[i]-1);
	}
	}

	if(screen) {
	fprintf(screen, "Printing symmetry functions...\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "Printing symmetry functions...\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	for(i=0; i<ne; i++) {
	if(screen) {
	fprintf(screen, "short range atomic symmetry functions element %2s\n", e[i]->name);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "short range atomic symmetry functions element %2s\n", e[i]->name);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	for(j=0; j<nsym[i]; j++) {
	G[i][j]->num = j;
	G[i][j]->info();
	}
	if(screen) fprintf(screen, "-----------------------------------------------------------------------\n");
	if(logfile) fprintf(logfile, "-----------------------------------------------------------------------\n");
	}

	return;

	}

	// read scaling.data and extract minimum, maximum and average SF values
	void PairRuNNer::read_scaling() {

	FILE *fp;
	char fname[LLINE+13];
	char line[LLINE];
	char *cerr;
	int ierr;
	int tmpi, tmpj;

	sprintf(fname, "%sscaling.data", dirRuNNer);
	fp = fopen(fname, "r");
	if(screen) fprintf(screen, "Reading scaling.data...\n");
	if(logfile) fprintf(logfile, "Reading scaling.data...\n");
	if(fp == NULL) {
	fprintf(stderr, "ERROR: Could not open scaling.data.\n");
	exit(EXIT_FAILURE);
	}
	else {
	for(int i=0; i<ne; i++) {
	for(int j=0; j<nsym[i]; j++) {
	cerr = fgets(line, LLINE, fp);
	ierr = sscanf(line, "%d %d %[^\n]", &tmpi, &tmpj, line);
	tmpi--;
	tmpj--;
	// save the scaling factor for later usage
	if(tmpi==i && tmpj==j) {
	ierr = sscanf(line, "%lf%lf%lf", &G[i][j]->min, &G[i][j]->max, &G[i][j]->ave);
	if (scalesym == 1 && centersym == 0) G[i][j]->sf = (scalemax - scalemin) / (G[i][j]->max - G[i][j]->min);
	else if(scalesym == 1 && centersym == 1) G[i][j]->sf = 1.0 / (G[i][j]->max - G[i][j]->min);
	else G[i][j]->sf = 1.0;
	}
	else {
	fprintf(stderr, "ERROR: scaling.data inconsistent with input.nn.\n");
	exit(EXIT_FAILURE);
	}
	}
	}
	fclose(fp);
	}
	if(screen) {
	fprintf(screen, "Finished reading scaling.data.\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "Finished reading scaling.data.\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}

	return;

	}

	// create symmetry function groups
	void PairRuNNer::create_symfunc_groups() {

	int kstart = 0;
	int** nsympg = NULL;
	int* tmptype = NULL;
	int* tmpe1 = NULL;
	int* tmpe2 = NULL;
	double* tmprc = NULL;
	double* tmpeta = NULL;
	double* tmpzeta = NULL;
	double* tmplambda = NULL;
	double* tmprs = NULL;

	if(comm->me == 0) {
	if(screen) {
	fprintf(screen, "Grouping symmetry functions...\n");
	}
	if(logfile) {
	fprintf(logfile, "Grouping symmetry functions...\n");
	}
	}

	nsymg = new int[ne];
	nsympg = new int*[ne];
	GG = new RuNNer_symfuncGroup**[ne];
	for(int i=0; i<ne; i++) {
	nsymg[i] = 0;
	nsympg[i] = new int[nsym[i]];
	for(int j=0; j<nsym[i]; j++) {
	nsympg[i][j] = 0;
	}
	}

	// detect number of symmetry function groups (assume symmetry functions are already sorted!)
	for(int i=0; i<ne; i++) {

	tmptype = new int[nsym[i]];
	tmprc = new double[nsym[i]];
	tmpeta = new double[nsym[i]];
	tmpzeta = new double[nsym[i]];
	tmplambda = new double[nsym[i]];
	tmprs = new double[nsym[i]];
	tmpe1 = new int[nsym[i]];
	tmpe2 = new int[nsym[i]];
	for(int j=0; j<nsym[i]; j++) {
	tmptype[j] = 0;
	tmprc[j] = 0.0;
	tmpeta[j] = 0.0;
	tmpzeta[j] = 0.0;
	tmplambda[j] = 0.0;
	tmprs[j] = 0.0;
	tmpe1[j] = 0;
	tmpe2[j] = 0;
	}


	for(int j=0; j<nsym[i]; j++) {
	if(G[i][j]->type == 2) {
	if(j == 0) {
	tmptype[0] = G[i][j]->type;
	tmprc [0] = G[i][j]->rc;
	tmpe1 [0] = G[i][j]->e1;
	nsymg[i]++;
	nsympg[i][0]++;
	continue;
	}
	for(int k=0; k<nsymg[i]; k++) {
	if((G[i][j]->type == tmptype[k]) &&
	(G[i][j]->rc == tmprc [k]) &&
	(G[i][j]->e1 == tmpe1 [k]) ) {
	nsympg[i][k]++;
	break;
	}
	if(k == nsymg[i] - 1) {
	tmptype[nsymg[i]] = G[i][j]->type;
	tmprc [nsymg[i]] = G[i][j]->rc ;
	tmpe1 [nsymg[i]] = G[i][j]->e1 ;
	nsympg[i][nsymg[i]]++;
	nsymg[i]++;
	break;
	}
	}
	}
	else if(G[i][j]->type == 3) {
	if(tmptype[nsymg[i] - 1] == 2) {
	kstart = nsymg[i];
	tmptype[kstart] = G[i][j]->type;
	tmprc [kstart] = G[i][j]->rc;
	tmpe1 [kstart] = G[i][j]->e1;
	tmpe2 [kstart] = G[i][j]->e2;
	nsymg[i]++;
	nsympg[i][kstart]++;
	continue;
	}
	for(int k=kstart; k<nsymg[i]; k++) {
	if((G[i][j]->type == tmptype[k]) &&
	(G[i][j]->rc == tmprc [k]) &&
	(G[i][j]->e1 == tmpe1 [k]) &&
	(G[i][j]->e2 == tmpe2 [k]) ) {
	nsympg[i][k]++;
	break;
	}
	if(k == nsymg[i] - 1) {
	tmptype[nsymg[i]] = G[i][j]->type;
	tmprc [nsymg[i]] = G[i][j]->rc ;
	tmpe1 [nsymg[i]] = G[i][j]->e1 ;
	tmpe2 [nsymg[i]] = G[i][j]->e2 ;
	nsympg[i][nsymg[i]]++;
	nsymg[i]++;
	break;
	}
	}
	}
	}

	// allocate symmetry function group objects
	GG[i] = new RuNNer_symfuncGroup*[nsymg[i]];
	for(int j=0; j<nsymg[i]; j++) GG[i][j] = new RuNNer_symfuncGroup(*this, j, nsympg[i][j]);

	// loop over symmetry functions again and setup group objects
	for(int j=0; j<nsym[i]; j++) {
	tmptype[j] = 0;
	tmprc[j] = 0.0;
	tmpeta[j] = 0.0;
	tmpzeta[j] = 0.0;
	tmplambda[j] = 0.0;
	tmprs[j] = 0.0;
	tmpe1[j] = 0;
	tmpe2[j] = 0;
	}

	nsymg[i] = 0;
	delete[] nsympg[i];
	nsympg[i] = new int[nsym[i]];
	for(int j=0; j<nsym[i]; j++) {
	nsympg[i][j] = 0;
	}


	for(int j=0; j<nsym[i]; j++) {
	if(G[i][j]->type == 2) {
	if(j == 0) {
	tmptype[0] = G[i][j]->type;
	tmprc [0] = G[i][j]->rc;
	tmpe1 [0] = G[i][j]->e1;
	GG[i][0]->setup(tmptype[0], i, tmpe1[0], tmpe2[0], tmprc[0]);
	GG[i][0]->setsc(0, G[i][j]->num, G[i][j]->min, G[i][j]->max, G[i][j]->ave, G[i][j]->sf);
	GG[i][0]->settype2(0, G[i][j]->eta, G[i][j]->rs);
	nsymg[i]++;
	nsympg[i][0]++;
	continue;
	}
	for(int k=0; k<nsymg[i]; k++) {
	if((G[i][j]->type == tmptype[k]) &&
	(G[i][j]->rc == tmprc [k]) &&
	(G[i][j]->e1 == tmpe1 [k]) ) {
	GG[i][k]->setsc(nsympg[i][k], G[i][j]->num, G[i][j]->min, G[i][j]->max, G[i][j]->ave, G[i][j]->sf);
	GG[i][k]->settype2(nsympg[i][k], G[i][j]->eta, G[i][j]->rs);
	nsympg[i][k]++;
	break;
	}
	if(k == nsymg[i] - 1) {
	tmptype[nsymg[i]] = G[i][j]->type;
	tmprc [nsymg[i]] = G[i][j]->rc ;
	tmpe1 [nsymg[i]] = G[i][j]->e1 ;
	GG[i][nsymg[i]]->setup(tmptype[nsymg[i]], i, tmpe1[nsymg[i]], tmpe2[nsymg[i]], tmprc[nsymg[i]]);
	GG[i][nsymg[i]]->setsc(nsympg[i][nsymg[i]], G[i][j]->num, G[i][j]->min, G[i][j]->max, G[i][j]->ave, G[i][j]->sf);
	GG[i][nsymg[i]]->settype2(nsympg[i][nsymg[i]], G[i][j]->eta, G[i][j]->rs);
	nsympg[i][nsymg[i]]++;
	nsymg[i]++;
	break;
	}
	}
	}
	else if(G[i][j]->type == 3) {
	if(tmptype[nsymg[i] - 1] == 2) {
	kstart = nsymg[i];
	tmptype[kstart] = G[i][j]->type;
	tmprc [kstart] = G[i][j]->rc;
	tmpe1 [kstart] = G[i][j]->e1;
	tmpe2 [kstart] = G[i][j]->e2;
	GG[i][kstart]->setup(tmptype[kstart], i, tmpe1[kstart], tmpe2[kstart], tmprc[kstart]);
	GG[i][kstart]->setsc(nsympg[i][kstart], G[i][j]->num, G[i][j]->min, G[i][j]->max, G[i][j]->ave, G[i][j]->sf);
	GG[i][kstart]->settype3(nsympg[i][kstart], G[i][j]->eta, G[i][j]->zeta, G[i][j]->lambda, G[i][j]->optpow);
	nsymg[i]++;
	nsympg[i][kstart]++;
	continue;
	}
	for(int k=kstart; k<nsymg[i]; k++) {
	if((G[i][j]->type == tmptype[k]) &&
	(G[i][j]->rc == tmprc [k]) &&
	(G[i][j]->e1 == tmpe1 [k]) &&
	(G[i][j]->e2 == tmpe2 [k]) ) {
	GG[i][k]->setsc(nsympg[i][k], G[i][j]->num, G[i][j]->min, G[i][j]->max, G[i][j]->ave, G[i][j]->sf);
	GG[i][k]->settype3(nsympg[i][k], G[i][j]->eta, G[i][j]->zeta, G[i][j]->lambda, G[i][j]->optpow);
	nsympg[i][k]++;
	break;
	}
	if(k == nsymg[i] - 1) {
	tmptype[nsymg[i]] = G[i][j]->type;
	tmprc [nsymg[i]] = G[i][j]->rc ;
	tmpe1 [nsymg[i]] = G[i][j]->e1 ;
	tmpe2 [nsymg[i]] = G[i][j]->e2 ;
	GG[i][nsymg[i]]->setup(tmptype[nsymg[i]], i, tmpe1[nsymg[i]], tmpe2[nsymg[i]], tmprc[nsymg[i]]);
	GG[i][nsymg[i]]->setsc(nsympg[i][nsymg[i]], G[i][j]->num, G[i][j]->min, G[i][j]->max, G[i][j]->ave, G[i][j]->sf);
	GG[i][nsymg[i]]->settype3(nsympg[i][nsymg[i]], G[i][j]->eta, G[i][j]->zeta, G[i][j]->lambda, G[i][j]->optpow);
	nsympg[i][nsymg[i]]++;
	nsymg[i]++;
	break;
	}
	}
	}
	}

	for(int j=0; j<nsymg[i]; j++) GG[i][j]->postset();

	delete[] tmptype ;
	delete[] tmprc ;
	delete[] tmpeta ;
	delete[] tmpzeta ;
	delete[] tmplambda;
	delete[] tmprs ;
	delete[] tmpe1 ;
	delete[] tmpe2 ;

	}

	if(comm->me == 0) {
	if(screen) {
	fprintf(screen, "Printing symmetry function groups...\n");
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "Printing symmetry function groups...\n");
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}

	for(int i=0; i<ne; i++) {
	if(screen) {
	fprintf(screen, "short range atomic symmetry functions groups of element %2s\n", e[i]->name);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "short range atomic symmetry functions groups of element %2s\n", e[i]->name);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	for(int j=0; j<nsymg[i]; j++) {
	GG[i][j]->info();
	}
	if(screen) {
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}
	}


	for(int i=0; i<ne; i++) delete[] nsympg[i];
	delete[] nsympg;

	return;

	}

	// calculate the maximum global cutoff for all symmetry functions
	double PairRuNNer::rc_max() {

	double rcmax = 0.0;

	for(int i=0; i<ne; i++) {
	for(int j=0; j<nsym[i]; j++) {
	rcmax = MAX(rcmax, G[i][j]->rc);
	}
	}
	if(comm->me==0) {
	if(screen) {
	fprintf(screen, "rcmax = %f\n", rcmax);
	fprintf(screen, "-----------------------------------------------------------------------\n");
	}
	if(logfile) {
	fprintf(logfile, "rcmax = %f\n", rcmax);
	fprintf(logfile, "-----------------------------------------------------------------------\n");
	}
	}

	return rcmax;

	}

	// soft cutoff function
	inline double PairRuNNer::fc(double r, double rcinv) {

	if(cutoff_type == 1) {
	return 0.5 * (cos(M_PI * r * rcinv) + 1.0);
	}
	else if(cutoff_type == 2) {
	double temp = tanh(1.0 - r * rcinv);
	return temp * temp * temp;
	}

	return -1;

	}

	// soft cutoff function derivative
	inline double PairRuNNer::dfc(double r, double rcinv) {

	if(cutoff_type == 1) {
	return -0.5 * M_PI * sin(M_PI * r * rcinv) * rcinv;
	}
	else if(cutoff_type == 2) {
	double temp = tanh(1.0 - r * rcinv);
	temp *= temp;
	return 3.0 * temp * (temp - 1.0) * rcinv;
	}

	return -1;

	}

	// optimized power function for integers (gsl_pow_int)
	inline double PairRuNNer::powint(double x, unsigned int n) {

	double value = 1.0;

	/* repeated squaring method
	* returns 0.0^0 = 1.0, so continuous in x
	*/
	do {
	if(n & 1) value = x; / for n odd */
	n >>= 1;
	x *= x;
	} while(n);

	return value;
	}


	// calculate one symmetry function
	double PairRuNNer::calc_G_dG(int ca, RuNNer_forces F, RuNNer_symfunc G, RuNNer_neighbor *N) {

	int num = G->num ;
	int type = G->type;
	//int ec = G->ec ;
	int e1 = G->e1 ;
	int e2 = G->e2 ;
	int optpow = G->optpow;
	double value = 0.0 ;
	double svalue = 0.0 ;
	double rc = G->rc ;
	double eta = G->eta ;
	double zeta = G->zeta ;
	double lambda = G->lambda;
	double rs = G->rs ;
	double min = G->min ;
	double max = G->max ;
	double ave = G->ave ;
	double sf = G->sf ;

	int Nej, Nek;
	int Nnum;
	double rcinv = 1.0 / rc;
	double rc2 = rc * rc;
	double rij, rik, rjk, rinvijik;
	double r2ij, r2ik;
	double dxij, dyij, dzij;
	double dxik, dyik, dzik;
	double dxjk, dyjk, dzjk;
	double costijk;
	double pexp, plambda, pfc, pdfc, pnorm, pzl;
	double pfcij, pfcik, pfcjk;
	double pdfcij;
	double p1, p2, p3;
	double p2etaplambda;
	double fg;
	double xptr, yptr, *zptr;

	// pointer to force class
	xptr = F->dGdx[num];
	yptr = F->dGdy[num];
	zptr = F->dGdz[num];

	// type 2 symmetry function
	if(type==2) {
	Nnum = N->num;
	for(int j=0; j<Nnum; j++) {
	Nej = N->e[j];
	rij = N->d[j];
	if(rij<=rc && e1 == Nej) {
	// energy calculation
	dxij = N->dx[j];
	dyij = N->dy[j];
	dzij = N->dz[j];
	pexp = exp(-eta * (rij - rs) * (rij - rs));
	pfc = fc(rij, rcinv);
	value += pexp * pfc;
	// force calculation
	pdfc = dfc(rij, rcinv);
	p1 = (pdfc - 2.0 * eta * (rij - rs) * pfc) * pexp / rij;
	dxij = p1 * dxij;
	dyij = p1 * dyij;
	dzij = p1 * dzij;
	// save force contributions in force class
	xptr[Nnum] += dxij;
	yptr[Nnum] += dyij;
	zptr[Nnum] += dzij;
	xptr[j] -= dxij;
	yptr[j] -= dyij;
	zptr[j] -= dzij;
	}
	}
	}

	// type 3 symmetry function
	if(type==3) {
	Nnum = N->num;
	pnorm = pow(2.0, 1.0 - zeta);
	pzl = zeta * lambda;
	for(int j=0; j<Nnum-1; j++) {
	Nej = N->e[j];
	rij = N->d[j];
	if( rij<=rc && (e1==Nej \|\| e2==Nej) ) {
	r2ij = rij * rij;
	pfcij = fc(rij, rcinv);
	pdfcij = dfc(rij, rcinv);
	for(int k=j+1; k<Nnum; k++) {
	Nek = N->e[k];
	if( (e1==Nej && e2==Nek) \|\|
	(e2==Nej && e1==Nek) ) {
	// energy calculation
	rik = N->d[k];
	if(rik <= rc) {
	dxjk = N->dx[k] - N->dx[j];
	dyjk = N->dy[k] - N->dy[j];
	dzjk = N->dz[k] - N->dz[j];
	rjk = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk;
	if(rjk <= rc2) {
	rjk = sqrt(rjk);
	rinvijik = 1.0 / rij / rik;
	r2ik = rik * rik;
	dxij = N->dx[j];
	dyij = N->dy[j];
	dzij = N->dz[j];
	dxik = N->dx[k];
	dyik = N->dy[k];
	dzik = N->dz[k];
	pfcik = fc(rik, rcinv);
	pfcjk = fc(rjk, rcinv);
	costijk = dxij * dxik + dyij * dyik + dzij * dzik;
	costijk *= rinvijik;
	pfc = pfcij * pfcik * pfcjk;
	pexp = exp(-eta * (r2ij + r2ik + rjk * rjk));
	plambda = 1.0 + lambda * costijk;
	if(plambda <= 0.0) fg = 0.0;
	else {
	if(optpow) fg = powint(plambda, optpow - 1) * pexp;
	else fg = pow(plambda, zeta - 1.0) * pexp;
	}
	value += fg * plambda * pfc;
	// force calculation
	fg *= pnorm;
	rinvijik *= pzl;
	costijk *= pzl;
	p2etaplambda = 2.0 * eta * plambda;
	p1 = pfc * (rinvijik - costijk / r2ij - p2etaplambda) + pfcik * pfcjk * pdfcij * plambda / rij;
	p2 = pfc * (rinvijik - costijk / r2ik - p2etaplambda) + pfcij * pfcjk * dfc(rik, rcinv) * plambda / rik;
	p3 = pfc * (rinvijik + p2etaplambda) - pfcij * pfcik * dfc(rjk, rcinv) * plambda / rjk;
	p1 *= fg;
	p2 *= fg;
	p3 *= fg;
	dxij *= p1;
	dyij *= p1;
	dzij *= p1;
	dxik *= p2;
	dyik *= p2;
	dzik *= p2;
	dxjk *= p3;
	dyjk *= p3;
	dzjk *= p3;
	// save force contributions in force class
	xptr[Nnum] += dxij + dxik;
	yptr[Nnum] += dyij + dyik;
	zptr[Nnum] += dzij + dzik;
	xptr[j] -= dxij + dxjk;
	yptr[j] -= dyij + dyjk;
	zptr[j] -= dzij + dzjk;
	xptr[k] -= dxik - dxjk;
	yptr[k] -= dyik - dyjk;
	zptr[k] -= dzik - dzjk;
	} // rjk <= rc
	} // rik <= rc
	} // elem
	} // k
	} // rij <= rc
	} // j
	value *= pnorm;
	} // type 3

	// scale and/or center symmetry functions
	if (scalesym == 1 && centersym == 0) svalue = (value - min) * sf + scalemin;
	else if(scalesym == 0 && centersym == 1) svalue = value - ave;
	else if(scalesym == 1 && centersym == 1) svalue = (value - ave) / (max - min);
	else svalue = value;

	// if SF value is out of boundaries give a warning
	if(value < min \|\| value > max) {
	if(showew == 1) {
	if(screen) {
	fprintf(screen, "PROC %d: ### EXTRAPOLATION WARNING ### TS: %ld ATOM: %d SYMFUNC: %d TYPE: %d VALUE: %f MIN: %f MAX: %f\n",
	comm->me, update->ntimestep, ca, num, type, value, min, max);
	}
	if(logfile) {
	fprintf(logfile, "PROC %d: ### EXTRAPOLATION WARNING ### TS: %ld ATOM: %d SYMFUNC: %d TYPE: %d VALUE: %f MIN: %f MAX: %f\n",
	comm->me, update->ntimestep, ca, num, type, value, min, max);
	}
	}
	//#ifdef _OPENMP
	// #pragma omp atomic
	ewcount++;
	ewcountsum++;
	//#endif
	}

	return svalue;

	}

	// calculate one symmetry function group
	void PairRuNNer::calc_G_dG_group(int ca, RuNNer_forces F, RuNNer_symfuncGroup GG, RuNNer_neighbor N, double& l_Gsv) {

	int numSF = GG->numSF ;
	int type = GG->type ;
	int e1 = GG->e1 ;
	int e2 = GG->e2 ;
	int*& index = GG->index ;
	double rc = GG->rc ;
	double*& min = GG->min ;
	double*& max = GG->max ;
	double*& ave = GG->ave ;
	double*& sf = GG->sf ;

	int*& optpow = GG->optpow;
	int*& etaind = GG->etaind;
	double*& eta = GG->eta ;
	double*& zeta = GG->zeta ;
	double*& lambda = GG->lambda;
	double*& rs = GG->rs ;

	int Nej, Nek;
	int Nnum;
	double rcinv = 1.0 / rc;
	double rc2 = rc * rc;
	double rij, rik, rjk, rinvijik;
	double r2ij, r2ik;
	double pr1, pr2, pr3;
	double dxij, dyij, dzij;
	double dxik, dyik, dzik;
	double dxjk, dyjk, dzjk;
	double costijk;
	double pexp, plambda, pfc, pdfc;
	double pfcij, pfcik, pfcjk;
	double pdfcij, pdfcik, pdfcjk;
	double p1, p2, p3;
	double p2etaplambda;
	double fg;

	double* value = NULL;
	double* svalue = NULL;
	double* pnorm = NULL;
	double r2sum = 0.0;
	double zl = 0.0;
	double vexp = 0.0;

	double* xptr = NULL;
	double* yptr = NULL;
	double* zptr = NULL;

	value = new double[numSF];
	svalue = new double[numSF];
	for(int i=0; i<numSF; i++) {
	value[i] = 0.0;
	svalue[i] = 0.0;
	}

	// type 2 symmetry function
	if(type == 2) {
	Nnum = N->num;
	for(int j=0; j<Nnum; j++) {
	Nej = N->e[j];
	rij = N->d[j];
	if(rij<=rc && e1 == Nej) {
	dxij = N->dx[j];
	dyij = N->dy[j];
	dzij = N->dz[j];
	pfc = fc(rij, rcinv);
	pdfc = dfc(rij, rcinv);
	for(int isym=0; isym<numSF; isym++) {
	// energy calculation
	pexp = exp(-eta[isym] * (rij - rs[isym]) * (rij - rs[isym]));
	value[isym] += pexp * pfc;
	// force calculation
	p1 = (pdfc - 2.0 * eta[isym] * (rij - rs[isym]) * pfc) * pexp / rij;
	// save force contributions in force class
	xptr = F->dGdx[index[isym]];
	yptr = F->dGdy[index[isym]];
	zptr = F->dGdz[index[isym]];
	xptr[Nnum] += p1 * dxij;
	yptr[Nnum] += p1 * dyij;
	zptr[Nnum] += p1 * dzij;
	xptr[j] -= p1 * dxij;
	yptr[j] -= p1 * dyij;
	zptr[j] -= p1 * dzij;
	}
	}
	}
	}

	// type 3 symmetry function
	if(type == 3) {
	pnorm = new double[numSF];
	for(int isym=0; isym<numSF; isym++) {
	pnorm[isym] = pow(2.0, 1.0 - zeta[isym]);
	}
	Nnum = N->num;
	for(int j=0; j<Nnum-1; j++) {
	Nej = N->e[j];
	rij = N->d[j];
	if( rij<=rc && (e1==Nej \|\| e2==Nej) ) {
	r2ij = rij * rij;
	pfcij = fc(rij, rcinv);
	pdfcij = dfc(rij, rcinv);
	dxij = N->dx[j];
	dyij = N->dy[j];
	dzij = N->dz[j];
	for(int k=j+1; k<Nnum; k++) {
	Nek = N->e[k];
	if( (e1==Nej && e2==Nek) \|\|
	(e2==Nej && e1==Nek) ) {
	// energy calculation
	rik = N->d[k];
	if(rik <= rc) {
	dxjk = N->dx[k] - dxij;
	dyjk = N->dy[k] - dyij;
	dzjk = N->dz[k] - dzij;
	rjk = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk;
	if(rjk <= rc2) {
	rjk = sqrt(rjk);
	rinvijik = 1.0 / (rij * rik);
	r2ik = rik * rik;
	dxik = N->dx[k];
	dyik = N->dy[k];
	dzik = N->dz[k];
	pfcik = fc(rik, rcinv);
	pfcjk = fc(rjk, rcinv);
	costijk = dxij * dxik + dyij * dyik + dzij * dzik;
	costijk *= rinvijik;
	pfc = pfcij * pfcik * pfcjk;
	pdfcik = dfc(rik, rcinv);
	pdfcjk = dfc(rjk, rcinv);
	r2sum = r2ij + r2ik + rjk * rjk;
	pr1 = pfcik * pfcjk * pdfcij / rij;
	pr2 = pfcij * pfcjk * pdfcik / rik;
	pr3 = pfcij * pfcik * pdfcjk / rjk;
	for(int isym=0; isym<numSF; isym++) {
	plambda = 1.0 + lambda[isym] * costijk;
	if(etaind[isym] == isym) {
	vexp = exp(-eta[isym] * r2sum);
	}
	if(plambda <= 0.0) fg = 0.0;
	else {
	if(optpow[isym]) fg = powint(plambda, optpow[isym] - 1) * vexp;
	else fg = pow(plambda, zeta[isym] - 1.0) * vexp;
	}
	value[isym] += fg * plambda * pfc;
	// force calculation
	fg *= pnorm[isym];
	zl = zeta[isym] * lambda[isym];
	p2etaplambda = 2.0 * eta[isym] * plambda / zl;
	p1 = fg * (pfc * zl * (rinvijik - costijk / r2ij - p2etaplambda) + pr1 * plambda);
	p2 = fg * (pfc * zl * (rinvijik - costijk / r2ik - p2etaplambda) + pr2 * plambda);
	p3 = fg * (pfc * zl * (rinvijik + p2etaplambda) - pr3 * plambda);
	// save force contributions in force class
	xptr = F->dGdx[index[isym]];
	yptr = F->dGdy[index[isym]];
	zptr = F->dGdz[index[isym]];
	xptr[Nnum] += p1 * dxij + p2 * dxik;
	yptr[Nnum] += p1 * dyij + p2 * dyik;
	zptr[Nnum] += p1 * dzij + p2 * dzik;
	xptr[j] -= p1 * dxij + p3 * dxjk;
	yptr[j] -= p1 * dyij + p3 * dyjk;
	zptr[j] -= p1 * dzij + p3 * dzjk;
	xptr[k] -= p2 * dxik - p3 * dxjk;
	yptr[k] -= p2 * dyik - p3 * dyjk;
	zptr[k] -= p2 * dzik - p3 * dzjk;
	} // isym
	} // rjk <= rc
	} // rik <= rc
	} // elem
	} // k
	} // rij <= rc
	} // j
	for(int isym=0; isym<numSF; isym++) {
	value[isym] *= pnorm[isym];
	}
	delete[] pnorm;
	} // type 3

	for(int isym=0; isym<numSF; isym++) {
	// scale and/or center symmetry functions
	if (scalesym == 1 && centersym == 0) svalue[isym] = (value[isym] - min[isym]) * sf[isym] + scalemin;
	else if(scalesym == 0 && centersym == 1) svalue[isym] = value[isym] - ave[isym];
	else if(scalesym == 1 && centersym == 1) svalue[isym] = (value[isym] - ave[isym]) / (max[isym] - min[isym]);
	else svalue[isym] = value[isym];

	// if SF value is out of boundaries give a warning
	if(value[isym] < min[isym] \|\| value[isym] > max[isym]) {
	if(showew == 1) {
	if(screen) {
	fprintf(screen, "PROC %d: ### EXTRAPOLATION WARNING ### TS: %ld ATOM: %d SYMFUNC: %d TYPE: %d VALUE: %f MIN: %f MAX: %f\n",
	comm->me, update->ntimestep, ca, index[isym], type, value[isym], min[isym], max[isym]);
	}
	if(logfile) {
	fprintf(logfile, "PROC %d: ### EXTRAPOLATION WARNING ### TS: %ld ATOM: %d SYMFUNC: %d TYPE: %d VALUE: %f MIN: %f MAX: %f\n",
	comm->me, update->ntimestep, ca, index[isym], type, value[isym], min[isym], max[isym]);
	}
	}
	ewcount++;
	ewcountsum++;
	}
	l_Gsv[index[isym]] = svalue[isym];
	}

	delete[] value;
	delete[] svalue;

	return;

	}

	// calculate all symmetry functions for all atoms
	void PairRuNNer::calc_all_symfunc() {

	int ca = 0;
	int ica = 0;
	int isym = 0;

	// loop over all elements
	for(int i=0; i<ne; i++) {
	// loop over all atoms of element i
	//#ifdef _OPENMP
	// #pragma omp parallel for private(ica, ca, isym)
	//#endif
	for(ica=0; ica<nlape[i]; ica++) {
	ca = laipe[i][ica];
	// calculate all SF of atom ca
	for(isym=0; isym<nsym[i]; isym++) {
	G_svalue[i][ica][isym] = calc_G_dG(la[ca], F[ca], G[i][isym], N[ca]);
	}
	}
	}

	if(ewcount > maxew) {
	if(screen) fflush(screen);
	if(logfile) fflush(logfile);
	fflush(stderr);
	usleep(1000000);
	fprintf(stderr, "ERROR: Too many extrapolation warnings.\n");
	exit(EXIT_FAILURE);
	}

	return;

	}

	// calculate all symmetry function groups for all atoms
	void PairRuNNer::calc_all_symfunc_group() {

	int ca = 0;
	int ica = 0;
	int isymg = 0;

	// loop over all elements
	for(int i=0; i<ne; i++) {
	// loop over all atoms of element i
	for(ica=0; ica<nlape[i]; ica++) {
	ca = laipe[i][ica];
	// calculate all SF groups of atom ca
	for(isymg=0; isymg<nsymg[i]; isymg++) {
	calc_G_dG_group(la[ca], F[ca], GG[i][isymg], N[ca], G_svalue[i][ica]);
	}
	}
	}

	if(ewcount > maxew) {
	if(screen) fflush(screen);
	if(logfile) fflush(logfile);
	fflush(stderr);
	usleep(1000000);
	fprintf(stderr, "ERROR: Too many extrapolation warnings.\n");
	exit(EXIT_FAILURE);
	}

	return;

	}

	// send neural network parameters to all processors
	void PairRuNNer::mpidistribute_nn() {

	// ALL
	MPI_Bcast(&ne, 1, MPI_INT, 0, world);
	MPI_Bcast(&nl, 1, MPI_INT, 0, world);
	MPI_Bcast(&scalesym, 1, MPI_INT, 0, world);
	MPI_Bcast(&centersym, 1, MPI_INT, 0, world);
	MPI_Bcast(&normnodes, 1, MPI_INT, 0, world);
	MPI_Bcast(&cutoff_type, 1, MPI_INT, 0, world);
	MPI_Bcast(&scalemin, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&scalemax, 1, MPI_DOUBLE, 0, world);

	// PROC >0
	if(comm->me>0) {
	e = new RuNNer_element*[ne];
	nsym = new int[ne];
	for(int i=0; i<ne; i++) {
	e[i] = new RuNNer_element(*this);
	nsym[i] = 0;
	}
	numnode = new int*[ne];
	actfunc = new int*[ne];
	for(int i=0; i<ne; i++) {
	numnode[i] = new int[nl+2];
	actfunc[i] = new int[nl+2];
	for(int j=0; j<nl+2; j++) {
	numnode[i][j] = 0;
	actfunc[i][j] = 0;
	}
	numnode[i][nl+1] = 1;
	actfunc[i][0] = 0;
	}
	}

	// ALL
	for(int i=0; i<ne; i++) {
	MPI_Bcast(&e[i]->nucch , 1, MPI_INT , 0, world);
	MPI_Bcast(&e[i]->atomic_energy, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&e[i]->name , 3, MPI_CHAR , 0, world);
	}
	MPI_Bcast(nsym, ne, MPI_INT, 0, world);
	for(int i=0; i<ne; i++) {
	MPI_Bcast(numnode[i], nl + 2, MPI_INT, 0, world);
	MPI_Bcast(actfunc[i], nl + 2, MPI_INT, 0, world);
	}

	return;

	}

	// send symmetry function parameters to all processors
	void PairRuNNer::mpidistribute_symfunc() {

	// PROC >0
	if(comm->me>0) {
	G = new RuNNer_symfunc**[ne];
	for(int i=0; i<ne; i++) {
	G[i] = new RuNNer_symfunc*[nsym[i]];
	for(int j=0; j<nsym[i]; j++) G[i][j] = new RuNNer_symfunc(*this);
	}
	}

	// ALL
	for(int i=0; i<ne; i++) {
	for(int j=0; j<nsym[i]; j++) {
	MPI_Bcast(&G[i][j]->num, 1, MPI_INT, 0, world);
	MPI_Bcast(&G[i][j]->type, 1, MPI_INT, 0, world);
	MPI_Bcast(&G[i][j]->ec, 1, MPI_INT, 0, world);
	MPI_Bcast(&G[i][j]->e1, 1, MPI_INT, 0, world);
	MPI_Bcast(&G[i][j]->e2, 1, MPI_INT, 0, world);
	MPI_Bcast(&G[i][j]->optpow, 1, MPI_INT, 0, world);
	MPI_Bcast(&G[i][j]->value, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->svalue, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->rc, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->eta, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->zeta, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->lambda, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->rs, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->min, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->max, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->ave, 1, MPI_DOUBLE, 0, world);
	MPI_Bcast(&G[i][j]->sf, 1, MPI_DOUBLE, 0, world);
	}
	}

	return;

	}

	// send neural network weights to all processors
	void PairRuNNer::mpidistribute_nnweights() {

	// ALL
	for(int i=0; i<ne; i++) {
	for(int j=0; j<nl+1; j++) {
	for(int k=0; k<L[i][j]->nprevnode; k++) {
	MPI_Bcast(L[i][j]->weight[k], L[i][j]->nnode, MPI_DOUBLE, 0, world);
	}
	MPI_Bcast(L[i][j]->bias, L[i][j]->nnode, MPI_DOUBLE, 0, world);
	}
	}

	return;

	}

	// sum up extrapolation warning counters of all processors
	void PairRuNNer::mpicollect_ew() {

	long ewcountsumglobal = 0;

	MPI_Allreduce(&ewcountsum, &ewcountsumglobal, 1, MPI_LONG, MPI_SUM, world);

	ewcountsum = ewcountsumglobal;

	return;

	}

	////////////////////////////////////////////////////////////////////////////////
	// PairRuNNer::RuNNer_atom member functions
	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::RuNNer_atom::RuNNer_atom(PairRuNNer& tmp): par(tmp) {

	e = 0;
	x = 0.0;
	y = 0.0;
	z = 0.0;
	en = 0.0;
	fx = 0.0;
	fy = 0.0;
	fz = 0.0;

	}

	PairRuNNer::RuNNer_atom::~RuNNer_atom() {}

	////////////////////////////////////////////////////////////////////////////////
	// PairRuNNer::RuNNer_element member functions
	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::RuNNer_element::RuNNer_element(PairRuNNer& x): par(x) {

	nucch = 0;
	strcpy(name, "??");
	atomic_energy = 0.0;

	}

	PairRuNNer::RuNNer_element::~RuNNer_element() {}

	PairRuNNer::RuNNer_element& PairRuNNer::RuNNer_element::operator=(const RuNNer_element& rhs) {

	if(this != &rhs) {
	nucch = rhs.nucch;
	strcpy(name, rhs.name);
	}

	return *this;

	}

	// get nuclear charge from name string
	int PairRuNNer::RuNNer_element::getnucch() {

	nucch = 0;

	if (strcmp(name, "H" ) == 0) nucch = 1;
	else if(strcmp(name, "He") == 0) nucch = 2;
	else if(strcmp(name, "Li") == 0) nucch = 3;
	else if(strcmp(name, "Be") == 0) nucch = 4;
	else if(strcmp(name, "B" ) == 0) nucch = 5;
	else if(strcmp(name, "C" ) == 0) nucch = 6;
	else if(strcmp(name, "N" ) == 0) nucch = 7;
	else if(strcmp(name, "O" ) == 0) nucch = 8;
	else if(strcmp(name, "F" ) == 0) nucch = 9;
	else if(strcmp(name, "Ne") == 0) nucch = 10;
	else if(strcmp(name, "Na") == 0) nucch = 11;
	else if(strcmp(name, "Mg") == 0) nucch = 12;
	else if(strcmp(name, "Al") == 0) nucch = 13;
	else if(strcmp(name, "Si") == 0) nucch = 14;
	else if(strcmp(name, "P" ) == 0) nucch = 15;
	else if(strcmp(name, "S" ) == 0) nucch = 16;
	else if(strcmp(name, "Cl") == 0) nucch = 17;
	else if(strcmp(name, "Ar") == 0) nucch = 18;
	else if(strcmp(name, "K" ) == 0) nucch = 19;
	else if(strcmp(name, "Ca") == 0) nucch = 20;
	else if(strcmp(name, "Sc") == 0) nucch = 21;
	else if(strcmp(name, "Ti") == 0) nucch = 22;
	else if(strcmp(name, "V" ) == 0) nucch = 23;
	else if(strcmp(name, "Cr") == 0) nucch = 24;
	else if(strcmp(name, "Mn") == 0) nucch = 25;
	else if(strcmp(name, "Fe") == 0) nucch = 26;
	else if(strcmp(name, "Co") == 0) nucch = 27;
	else if(strcmp(name, "Ni") == 0) nucch = 28;
	else if(strcmp(name, "Cu") == 0) nucch = 29;
	else if(strcmp(name, "Zn") == 0) nucch = 30;
	else if(strcmp(name, "Ga") == 0) nucch = 31;
	else if(strcmp(name, "Ge") == 0) nucch = 32;
	else if(strcmp(name, "As") == 0) nucch = 33;
	else if(strcmp(name, "Se") == 0) nucch = 34;
	else if(strcmp(name, "Br") == 0) nucch = 35;
	else if(strcmp(name, "Kr") == 0) nucch = 36;
	else if(strcmp(name, "Rb") == 0) nucch = 37;
	else if(strcmp(name, "Sr") == 0) nucch = 38;
	else if(strcmp(name, "Y" ) == 0) nucch = 39;
	else if(strcmp(name, "Zr") == 0) nucch = 40;
	else if(strcmp(name, "Nb") == 0) nucch = 41;
	else if(strcmp(name, "Mo") == 0) nucch = 42;
	else if(strcmp(name, "Tc") == 0) nucch = 43;
	else if(strcmp(name, "Ru") == 0) nucch = 44;
	else if(strcmp(name, "Rh") == 0) nucch = 45;
	else if(strcmp(name, "Pd") == 0) nucch = 46;
	else if(strcmp(name, "Ag") == 0) nucch = 47;
	else if(strcmp(name, "Cd") == 0) nucch = 48;
	else if(strcmp(name, "In") == 0) nucch = 49;
	else if(strcmp(name, "Sn") == 0) nucch = 50;
	else if(strcmp(name, "Sb") == 0) nucch = 51;
	else if(strcmp(name, "Te") == 0) nucch = 52;
	else if(strcmp(name, "I" ) == 0) nucch = 53;
	else if(strcmp(name, "Xe") == 0) nucch = 54;
	else if(strcmp(name, "Cs") == 0) nucch = 55;
	else if(strcmp(name, "Ba") == 0) nucch = 56;
	else if(strcmp(name, "La") == 0) nucch = 57;
	else if(strcmp(name, "Ce") == 0) nucch = 58;
	else if(strcmp(name, "Pr") == 0) nucch = 59;
	else if(strcmp(name, "Nd") == 0) nucch = 60;
	else if(strcmp(name, "Pm") == 0) nucch = 61;
	else if(strcmp(name, "Sm") == 0) nucch = 62;
	else if(strcmp(name, "Eu") == 0) nucch = 63;
	else if(strcmp(name, "Gd") == 0) nucch = 64;
	else if(strcmp(name, "Tb") == 0) nucch = 65;
	else if(strcmp(name, "Dy") == 0) nucch = 66;
	else if(strcmp(name, "Ho") == 0) nucch = 67;
	else if(strcmp(name, "Er") == 0) nucch = 68;
	else if(strcmp(name, "Tm") == 0) nucch = 69;
	else if(strcmp(name, "Yb") == 0) nucch = 70;
	else if(strcmp(name, "Lu") == 0) nucch = 71;
	else if(strcmp(name, "Hf") == 0) nucch = 72;
	else if(strcmp(name, "Ta") == 0) nucch = 73;
	else if(strcmp(name, "W" ) == 0) nucch = 74;
	else if(strcmp(name, "Re") == 0) nucch = 75;
	else if(strcmp(name, "Os") == 0) nucch = 76;
	else if(strcmp(name, "Ir") == 0) nucch = 77;
	else if(strcmp(name, "Pt") == 0) nucch = 78;
	else if(strcmp(name, "Au") == 0) nucch = 79;
	else if(strcmp(name, "Hg") == 0) nucch = 80;
	else if(strcmp(name, "Tl") == 0) nucch = 81;
	else if(strcmp(name, "Pb") == 0) nucch = 82;
	else if(strcmp(name, "Bi") == 0) nucch = 83;
	else if(strcmp(name, "Po") == 0) nucch = 84;
	else if(strcmp(name, "At") == 0) nucch = 85;
	else if(strcmp(name, "Rn") == 0) nucch = 86;
	else if(strcmp(name, "Fr") == 0) nucch = 87;
	else if(strcmp(name, "Ra") == 0) nucch = 88;
	else if(strcmp(name, "Ac") == 0) nucch = 89;
	else if(strcmp(name, "Th") == 0) nucch = 90;
	else if(strcmp(name, "Pa") == 0) nucch = 91;
	else if(strcmp(name, "U" ) == 0) nucch = 92;
	else if(strcmp(name, "Np") == 0) nucch = 93;
	else if(strcmp(name, "Pu") == 0) nucch = 94;
	else if(strcmp(name, "Am") == 0) nucch = 95;
	else if(strcmp(name, "Cm") == 0) nucch = 96;
	else if(strcmp(name, "Bk") == 0) nucch = 97;
	else if(strcmp(name, "Cf") == 0) nucch = 98;
	else if(strcmp(name, "Es") == 0) nucch = 99;
	else if(strcmp(name, "Fm") == 0) nucch = 100;
	else if(strcmp(name, "Md") == 0) nucch = 101;
	else if(strcmp(name, "No") == 0) nucch = 102;
	else {
	fprintf(stderr, "ERROR: Unknown element %s.\n", name);
	exit(EXIT_FAILURE);
	}

	return nucch;

	}

	////////////////////////////////////////////////////////////////////////////////
	// PairRuNNer::RuNNer_forces member functions
	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::RuNNer_forces::RuNNer_forces(PairRuNNer& x, int nlayer, int *nnlayer, int element, int ai, int nn): par(x) {

	e = element;
	nlf = nlayer + 1;
	nsym = nnlayer[0];
	indatom = ai;
	num = nn;
	nnode = new int[nlf];
	dEdG = new double[nsym];
	dfdx = new double*[nlf];
	for(int i=0; i<nlf; i++) {
	nnode[i] = nnlayer[i+1];
	dfdx[i] = new double[nnlayer[i+1]];
	for(int j=0; j<nnlayer[i+1]; j++) {
	dfdx[i][j] = 0.0;
	}
	}
	dGdx = new double*[nsym];
	dGdy = new double*[nsym];
	dGdz = new double*[nsym];
	for(int i=0; i<nsym; i++) {
	dEdG[i] = 0.0;
	dGdx[i] = new double[num+1];
	dGdy[i] = new double[num+1];
	dGdz[i] = new double[num+1];
	for(int j=0; j<num+1; j++) {
	dGdx[i][j] = 0.0;
	dGdy[i][j] = 0.0;
	dGdz[i][j] = 0.0;
	}
	}

	}

	PairRuNNer::RuNNer_forces::~RuNNer_forces() {

	delete[] nnode;
	for(int i=0; i<nlf; i++) {
	delete[] dfdx[i];
	}
	delete[] dfdx;
	for(int i=0; i<nsym; i++) {
	delete[] dGdx[i];
	delete[] dGdy[i];
	delete[] dGdz[i];
	}
	delete[] dGdx;
	delete[] dGdy;
	delete[] dGdz;
	delete[] dEdG;

	}

	// get saved df/dx from layer class
	void PairRuNNer::RuNNer_forces::get_dfdx(RuNNer_layer **L) {

	for(int i=0; i<nlf; i++) {
	for(int j=0; j<nnode[i]; j++) {
	dfdx[i][j] = L[i]->dfdx[j];
	}
	}

	return;

	}

	// calculate dE/dG from saved values in layer class
	void PairRuNNer::RuNNer_forces::calc_dEdG(RuNNer_layer **L) {

	double tmpinner, tmpouter;

	tmpinner = new double*[nlf-1];
	tmpouter = new double*[nlf-1];
	for(int i=0; i<nlf-1; i++) {
	tmpinner[i] = new double[nnode[i]];
	tmpouter[i] = new double[nnode[i+1]];
	}

	for(int k=0; k<nsym; k++) {
	for(int i=0; i<nnode[0]; i++) {
	tmpinner[0][i] = dfdx[0][i] * L[0]->weight[k][i];
	if(par.normnodes) tmpinner[0][i] /= L[0]->nprevnode;
	}
	for(int l=0; l<nlf-1; l++) {
	for(int i=0; i<nnode[l+1]; i++) {
	tmpouter[l][i] = 0.0;
	for(int j=0; j<nnode[l]; j++) {
	tmpouter[l][i] += L[l+1]->weight[j][i] * tmpinner[l][j];
	}
	tmpouter[l][i] *= dfdx[l+1][i];
	if(par.normnodes) tmpouter[l][i] /= L[l+1]->nprevnode;
	if(l<nlf-2) tmpinner[l+1][i] = tmpouter[l][i];
	}
	}
	dEdG[k] = tmpouter[nlf-2][0];
	}

	for(int i=0; i<nlf-1; i++) {
	delete[] tmpinner[i];
	delete[] tmpouter[i];
	}
	delete[] tmpinner;
	delete[] tmpouter;

	return;

	}

	////////////////////////////////////////////////////////////////////////////////
	// PairRuNNer::RuNNer_layer member functions
	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::RuNNer_layer::RuNNer_layer(PairRuNNer& x, int a, int b, int c): par(x) {

	nprevnode = a;
	nnode = b;
	af = c;
	input = new double[a];
	node = new double[b];
	dfdx = new double[b];
	bias = new double[b];
	weight = new double*[a];
	for(int i=0; i<a; i++) {
	input[i] = 0.0;
	weight[i] = new double[b];
	for(int j=0; j<b; j++) weight[i][j] = 0.0;
	}
	for(int i=0; i<b; i++) {
	node[i] = 0.0;
	dfdx[i] = 0.0;
	bias[i] = 0.0;
	}

	}

	PairRuNNer::RuNNer_layer::~RuNNer_layer() {

	delete[] input;
	delete[] node;
	delete[] dfdx;
	delete[] bias;
	for(int i=0; i<nprevnode; i++) {
	delete[] weight[i];
	}
	delete[] weight;

	}

	// get NN bias and weights from weights.???.data
	void PairRuNNer::RuNNer_layer::get_weights(FILE *fp) {

	char line[LLINE];
	int ierr;
	char *cerr;

	// read weights
	for(int w1=0; w1<nprevnode; w1++) {
	for(int w2=0; w2<nnode; w2++) {
	cerr = fgets(line, LLINE, fp);
	ierr = sscanf(line, "%lf", &weight[w1][w2]);
	}
	}
	// read bias
	for(int b=0; b<nnode; b++) {
	cerr = fgets(line, LLINE, fp);
	ierr = sscanf(line, "%lf", &bias[b]);
	}

	return;

	}

	// calc all nodes of layer
	void PairRuNNer::RuNNer_layer::calc_all_nodes() {

	for(int i=0; i<nnode; i++) {
	calc_one_node(i);
	}

	return;

	}

	// calc one node of layer
	void PairRuNNer::RuNNer_layer::calc_one_node(int inode) {

	double nv = 0.0;

	for(int j=0; j<nprevnode; j++) {
	nv += weight[j][inode] * input[j];
	}
	nv += bias[inode];
	if(par.normnodes) nv /= nprevnode;

	// save node values and derivatives for later usage
	if(af==1) {
	dfdx[inode] = 1.0;
	node[inode] = nv;
	}
	else if(af==2) {
	dfdx[inode] = 1.0 - pow(tanh(nv), 2.0);
	node[inode] = tanh(nv);
	}
	else {
	fprintf(stderr, "ERROR: Unknown activation function.\n");
	exit(EXIT_FAILURE);
	}

	return;

	}

	// get input node values from previous layer
	void PairRuNNer::RuNNer_layer::get_input(RuNNer_layer *lprev) {

	for(int i=0; i<nprevnode; i++) {
	input[i] = lprev->node[i];
	}

	return;

	}

	////////////////////////////////////////////////////////////////////////////////
	// PairRuNNer::RuNNer_neighbor member functions
	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::RuNNer_neighbor::RuNNer_neighbor(PairRuNNer& x): par(x) {

	iatom = 0;
	eatom = 0;
	num = 0;
	for(int i=0; i<MAXNEIGH; i++) {
	ilal[i] = 0;
	e[i] = 0;
	d[i] = 0.0;
	dx[i] = 0.0;
	dy[i] = 0.0;
	dz[i] = 0.0;
	}

	}

	PairRuNNer::RuNNer_neighbor::~RuNNer_neighbor() {}

	// calculate neighbor list using LAMMPS neighbor list
	int PairRuNNer::RuNNer_neighbor::calcl(RuNNer_atom **a, int ii) {

	int i = par.list->ilist[ii];
	int ila = par.la[i];
	int j = 0;
	double tmpd;
	double tmpdx, tmpdy, tmpdz;
	int numneigh = par.list->numneigh[i];
	int* firstneigh = par.list->firstneigh[i];

	eatom = a[ila]->e;

	for(int jj=0; jj<numneigh; jj++) {
	j = firstneigh[jj];
	tmpdx = a[ila]->x - CONST_A2B * par.atom->x[j][0];
	tmpdy = a[ila]->y - CONST_A2B * par.atom->x[j][1];
	tmpdz = a[ila]->z - CONST_A2B * par.atom->x[j][2];
	tmpd = sqrt(tmpdx * tmpdx + tmpdy * tmpdy + tmpdz * tmpdz);
	// if distance <= maximal cutoff add atom to neighbor list
	if(tmpd <= par.rcmax) {
	if(num >= MAXNEIGH) {
	fprintf(stderr, "ERROR: neighbor list array too small, increase MAXNEIGH and recompile.\n");
	exit(EXIT_FAILURE);
	}
	ilal[num] = j;
	//ilal2[j] = num;
	e[num] = par.transform_atom_type(par.atom->type[j])-1; //HELLSTROM
	d[num] = tmpd;
	dx[num] = tmpdx;
	dy[num] = tmpdy;
	dz[num] = tmpdz;
	num++;
	}
	}

	return num;

	}

	// calculate neighbor list
	int PairRuNNer::RuNNer_neighbor::calc(RuNNer_atom **a, int i) {

	int j;
	double tmpd;
	double tmpdx, tmpdy, tmpdz;

	eatom = a[i]->e;

	for(j=0; j<par.atom->nlocal+par.atom->nghost; j++) {
	if( !( i==par.atom->tag[j]-1 && j<par.nla ) ) {
	tmpdx = a[i]->x - CONST_A2B * par.atom->x[j][0];
	tmpdy = a[i]->y - CONST_A2B * par.atom->x[j][1];
	tmpdz = a[i]->z - CONST_A2B * par.atom->x[j][2];
	tmpd = sqrt(tmpdx * tmpdx + tmpdy * tmpdy + tmpdz * tmpdz);
	// if distance <= maximal cutoff add atom to neighbor list
	if(tmpd <= par.rcmax) {
	if(num >= MAXNEIGH) {
	fprintf(stderr, "ERROR: neighbor list array too small, increase MAXNEIGH and recompile.\n");
	exit(EXIT_FAILURE);
	}
	ilal[num] = j;
	//ilal2[j] = num;
	e[num] = par.transform_atom_type(par.atom->type[j])-1; //HELLSTROM
	d[num] = tmpd;
	dx[num] = tmpdx;
	dy[num] = tmpdy;
	dz[num] = tmpdz;
	num++;
	}
	}
	}

	return num;

	}

	// find LAMMPS atom list index in RuNNer neighbor list
	int PairRuNNer::RuNNer_neighbor::find_nindex(int n) {

	for(int i=0; i<num; i++) {
	if(ilal[i]==n) return i;
	}

	return -1;

	}

	////////////////////////////////////////////////////////////////////////////////
	// PairRuNNer::RuNNer_symfunc member functions
	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::RuNNer_symfunc::RuNNer_symfunc(PairRuNNer& x): par(x) {

	num = 0;
	type = 0;
	ec = 0;
	e1 = 0;
	e2 = 0;
	optpow = 0;
	value = 0.0;
	svalue = 0.0;
	rc = 0.0;
	eta = 0.0;
	zeta = 0.0;
	lambda = 0.0;
	rs = 0.0;
	min = 0.0;
	max = 0.0;
	ave = 0.0;
	sf = 0.0;

	}

	PairRuNNer::RuNNer_symfunc::~RuNNer_symfunc() {}

	PairRuNNer::RuNNer_symfunc& PairRuNNer::RuNNer_symfunc::operator=(const RuNNer_symfunc& rhs) {

	if(this != &rhs) {
	num = rhs.num;
	type = rhs.type;
	ec = rhs.ec;
	e1 = rhs.e1;
	e2 = rhs.e2;
	optpow = rhs.optpow;
	value = rhs.value;
	svalue = rhs.svalue;
	rc = rhs.rc;
	eta = rhs.eta;
	zeta = rhs.zeta;
	lambda = rhs.lambda;
	rs = rhs.rs;
	min = rhs.min;
	max = rhs.max;
	ave = rhs.ave;
	sf = rhs.sf;
	}

	return *this;

	}

	// set up one symmetry function
	void PairRuNNer::RuNNer_symfunc::setup(char *line) {

	int tmp;
	char tmpename[3];
	char tmpe1name[3];
	char tmpe2name[3];

	// analyze char string passed from analyze_input_nn()
	sscanf(line, "%s %d", tmpename, &type);
	ec = par.find_element(tmpename);
	if(type==2) {
	sscanf(line, "%s %s %s %lf %lf %lf", tmpe1name, &eta, &rs, &rc);
	e1 = par.find_element(tmpe1name);
	}
	else if(type==3) {
	sscanf(line, "%s %s %s %s %lf %lf %lf %lf", tmpe1name, tmpe2name, &eta, &lambda, &zeta, &rc);
	e1 = par.find_element(tmpe1name);
	e2 = par.find_element(tmpe2name);
	if(e1>e2) {
	tmp = e1;
	e1 = e2;
	e2 = tmp;
	}
	if(fabs(zeta - round(zeta)) < SMALL_DOUBLE) {
	optpow = round(zeta);
	}
	}

	return;

	}

	// info output
	void PairRuNNer::RuNNer_symfunc::info() {

	if(type==2) {
	if(par.screen) fprintf(par.screen, "%4d %2s %2d %2s %7.3f %7.3f %7.3f\n", num + 1, par.e[ec]->name, type, par.e[e1]->name, eta, rs, rc);
	if(par.logfile) fprintf(par.logfile, "%4d %2s %2d %2s %7.3f %7.3f %7.3f\n", num + 1, par.e[ec]->name, type, par.e[e1]->name, eta, rs, rc);
	}
	else if(type==3) {
	if(par.screen) fprintf(par.screen, "%4d %2s %2d %2s %2s %7.3f %7.3f %7.3f %7.3f\n", num + 1, par.e[ec]->name, type, par.e[e1]->name, par.e[e2]->name, eta, lambda, zeta, rc);
	if(par.logfile) fprintf(par.logfile, "%4d %2s %2d %2s %2s %7.3f %7.3f %7.3f %7.3f\n", num + 1, par.e[ec]->name, type, par.e[e1]->name, par.e[e2]->name, eta, lambda, zeta, rc);
	}

	return;

	}

	////////////////////////////////////////////////////////////////////////////////
	// PairRuNNer::RuNNer_symfuncGroup member functions
	////////////////////////////////////////////////////////////////////////////////

	PairRuNNer::RuNNer_symfuncGroup::RuNNer_symfuncGroup(PairRuNNer& x, int l_num, int l_numSF): par(x) {

	num = l_num;
	numSF = l_numSF;
	type = 0;
	ec = 0;
	e1 = 0;
	e2 = 0;
	index = NULL;
	optpow = NULL;
	etaind = NULL;
	rc = 0.0;
	eta = NULL;
	zeta = NULL;
	lambda = NULL;
	rs = NULL;
	min = NULL;
	max = NULL;
	ave = NULL;
	sf = NULL;

	}

	PairRuNNer::RuNNer_symfuncGroup::~RuNNer_symfuncGroup() {

	if(index != NULL) delete[] index ;
	if(optpow != NULL) delete[] optpow;
	if(etaind != NULL) delete[] etaind;
	if(eta != NULL) delete[] eta ;
	if(zeta != NULL) delete[] zeta ;
	if(lambda != NULL) delete[] lambda;
	if(rs != NULL) delete[] rs ;
	if(min != NULL) delete[] min ;
	if(max != NULL) delete[] max ;
	if(ave != NULL) delete[] ave ;
	if(sf != NULL) delete[] sf ;

	}

	void PairRuNNer::RuNNer_symfuncGroup::setup(int l_type, int l_ec, int l_e1, int l_e2, double l_rc) {

	// SF group parameters
	type = l_type;
	ec = l_ec;
	e1 = l_e1;
	e2 = l_e2;
	rc = l_rc;

	// allocate type-specific arrays
	if(type == 2) {
	eta = new double[numSF];
	rs = new double[numSF];
	for(int i=0; i<numSF; i++) {
	eta [i] = 0.0;
	rs [i] = 0.0;
	}
	}
	else if(type == 3) {
	optpow = new int [numSF];
	etaind = new int [numSF];
	eta = new double[numSF];
	zeta = new double[numSF];
	lambda = new double[numSF];
	for(int i=0; i<numSF; i++) {
	optpow[i] = 0 ;
	etaind[i] = 0 ;
	eta [i] = 0.0;
	zeta [i] = 0.0;
	lambda[i] = 0.0;
	}

	}

	// allocate general SF arrays
	index = new int [numSF];
	min = new double[numSF];
	max = new double[numSF];
	ave = new double[numSF];
	sf = new double[numSF];
	for(int i=0; i<numSF; i++) {
	index [i] = 0 ;
	min [i] = 0.0;
	max [i] = 0.0;
	ave [i] = 0.0;
	sf [i] = 0.0;
	}

	return;

	}

	void PairRuNNer::RuNNer_symfuncGroup::info() {

	if(type == 2) {
	if(par.screen) {
	fprintf(par.screen, "%4d %2s %2d %2s * * %7.3f %4d\n", num + 1, par.e[ec]->name, type, par.e[e1]->name, rc, numSF);
	for(int i=0; i<numSF; i++) {
	fprintf(par.screen, " %7.3f %7.3f %4d %4d\n", eta[i], rs[i], i + 1, index[i] + 1);
	}
	}
	if(par.logfile) {
	fprintf(par.logfile, "%4d %2s %2d %2s * * %7.3f %4d\n", num + 1, par.e[ec]->name, type, par.e[e1]->name, rc, numSF);
	for(int i=0; i<numSF; i++) {
	fprintf(par.logfile, " %7.3f %7.3f %4d %4d\n", eta[i], rs[i], i + 1, index[i] + 1);
	}
	}
	}
	else if(type == 3) {
	if(par.screen) {
	fprintf(par.screen, "%4d %2s %2d %2s %2s * * * %7.3f %4d\n", num + 1, par.e[ec]->name, type, par.e[e1]->name, par.e[e2]->name, rc, numSF);
	for(int i=0; i<numSF; i++) {
	fprintf(par.screen, " %7.3f %7.3f %7.3f %4d %4d %4d\n", eta[i], lambda[i], zeta[i], i + 1, index[i] + 1, etaind[i] + 1);
	}
	}
	if(par.logfile) {
	fprintf(par.logfile, "%4d %2s %2d %2s %2s * * * %7.3f %4d\n", num + 1, par.e[ec]->name, type, par.e[e1]->name, par.e[e2]->name, rc, numSF);
	for(int i=0; i<numSF; i++) {
	fprintf(par.logfile, " %7.3f %7.3f %7.3f %4d %4d %4d\n", eta[i], lambda[i], zeta[i], i + 1, index[i] + 1, etaind[i] + 1);
	}
	}
	}

	return;

	}

	void PairRuNNer::RuNNer_symfuncGroup::setsc(int i, int l_index, double l_min, double l_max, double l_ave, double l_sf) {

	index [i] = l_index ;
	min [i] = l_min ;
	max [i] = l_max ;
	ave [i] = l_ave ;
	sf [i] = l_sf ;

	return;

	}

	void PairRuNNer::RuNNer_symfuncGroup::settype2(int i, double l_eta, double l_rs) {

	eta [i] = l_eta ;
	rs [i] = l_rs ;

	return;

	}

	void PairRuNNer::RuNNer_symfuncGroup::settype3(int i, double l_eta, double l_zeta, double l_lambda, int l_optpow) {

	eta [i] = l_eta ;
	zeta [i] = l_zeta ;
	lambda[i] = l_lambda;
	optpow[i] = l_optpow;

	return;

	}

	void PairRuNNer::RuNNer_symfuncGroup::postset() {

	if(type == 3) {

	int index = 0;
	double etatmp = 0.0;

	index = 0;
	etaind[0] = index;
	etatmp = eta[0];
	for(int i=1; i<numSF; i++) {
	if(etatmp == eta[i]) {
	etaind[i] = index;
	}
	else {
	index = i;
	etaind[i] = index;
	etatmp = eta[i];
	}
	}
	}

	return;

	}

pair_runner.cppNo OneTemporaryActions

File Metadata

pair_runner.cppView Options

Event Timeline

pair_runner.cpp
No OneTemporary
Actions

pair_runner.cpp
View Options