Homec4science
Diffusion lammps (master) Browse

lammps/src/NEURONET/pair_neuronet.cpp
7a23f3894a40master

Last Change
Raw File

pair_neuronet.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
This file is written by Till Junge <till.junge@epfl.ch>
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_neuronet.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
//#include "respa.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include <fstream>
extern "C" {
void __nn_MOD_force_nn(int * namax, int * natm, int* tag, double* ra,
int * nnmax, double * force, double * vatom, double *h,
double *hi, double * tcom, int * nb, int * nbmax,
int * lsb, int * nex, int * lsrc, int * myparity, int * nn,
double * sv, double * rc, int *lspr, int * mpi_world,
int * myid, double * epi, double * epot, bool * vflag_atom);
void __nn_MOD_read_params(int * myid, int * mpi_world, double * rcin,
double * rc3);
double __nn_MOD_dsigmoid(double * x);
}
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairNeuroNet::PairNeuroNet(LAMMPS *lmp) : Pair(lmp)
{
respa_enable = 1;
writedata = 1;
this->ncnst_type[0]= 2; // Gaussian
this->ncnst_type[1]= 1; // cosine
this->ncnst_type[2]= 1; // polynomial
this->ncnst_type[3]= 2; // Morse
this->ncnst_type[100]= 1; // angular
for (int i = 0; i < 100; i++) {
this->ncomb_type[i]= 2; // pair
this->ncomb_type[100+i] = 3; // triplet
}
}
/* ---------------------------------------------------------------------- */
PairNeuroNet::~PairNeuroNet()
{
if (allocated) {
delete[] this->itype;
delete[] this->cnst;
delete[] this->iaddr2;
delete[] this->iaddr3;
delete[] this->wgt11;
delete[] this->wgt12;
delete[] this->wgt21;
delete[] this->wgt22;
delete[] this->wgt23;
}
}
/* ---------------------------------------------------------------------- */
void PairNeuroNet::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
int nmax = atom->nmax;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
inum = list->inum; //number of atoms that have a neighbor list
ilist = list->ilist; //indices of those atoms
numneigh = list->numneigh; // number of their neighbors
firstneigh = list->firstneigh; // pointer to first neighbor index (indices contiguous
// loop over neighbors of my atoms
int nnmax = 0;
for (int ii = 0; ii < inum; ++ii) {
i = ilist[ii];
nnmax = (nnmax > numneigh[i]) ? nnmax : numneigh[i];
}
// ryo neighlist is matrix
int fort_neighbor[(nnmax+1) * nmax] = {0};
int fort_type[nmax] = {0};
double fort_x [3*nmax] = {0};
double fort_f [3*nmax] = {0};
double fort_vatom [6*nmax] = {0};
double fort_eatom [nmax] = {0};
for (int ii = 0; ii < inum; ++ii) {
i = ilist[ii];
for (int dir = 0; dir < 3; ++dir) {
fort_x[dir+3*ii] = x[i][dir];
}
int *neigh_ptr = firstneigh[i];
int nb_neigh = numneigh[i];
fort_neighbor[0 + (nnmax+1)*ii] = nb_neigh;
for (int j = 0; j < nb_neigh; j++) {
fort_neighbor[j+1 + (nnmax+1)*ii] = neigh_ptr[j];
}
fort_type[ii] = type[i];
}
double h[] = {1,0,0,
0,1,0,
0,0,1};
double dummy = 0;
int idummy = 0;
double delta_eng_vdwl;
bool fort_vflag_atom = bool(this->vflag_atom);
__nn_MOD_force_nn(&nmax, &inum, fort_type, fort_x, &nnmax, fort_f,
fort_vatom, h, h, &dummy, &nghost, &idummy, &idummy, &idummy,
&idummy, &idummy, &idummy, &dummy, &rcin, fort_neighbor, &idummy,
&idummy, fort_eatom, &delta_eng_vdwl, &fort_vflag_atom);
this->eng_vdwl += delta_eng_vdwl;
// rewrite results into lammps form
for (int ii = 0; ii < inum; ++ii) {
i = ilist[ii];
eatom[i] += fort_eatom[ii];
for (int dir = 0; dir < 3; ++dir) {
f[i][dir] += fort_f[dir+3*ii];
}
for (int dir = 0; dir < 6; ++dir) {
vatom[i][dir] += fort_vatom[dir+6*ii];
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairNeuroNet::settings(int narg, char **arg) // read_params
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairNeuroNet::coeff(int narg, char **arg)
{
if (narg < 5)
error->all(FLERR,"Not ennough arguments, should be * * constants_input_file params_input_file ");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// Read constants input file
std::ifstream const_file(arg[2]);
//TODO: check that file exist etc
int dummy;
const_file >> this->nl >> this->nsp >> dummy;
// dummy contains number of input nodes
this->nhl.push_back(dummy);
for (int i = 0; i < nl; ++i) {
// read the number of nodes per hidden layer
const_file >> dummy;
this->nhl.push_back(dummy);
}
// last layer always has one node
nhl.push_back(1);
// allocate arrays for consts
this->itype = new int[this->nhl[0]];
this->cnst = new double[this->max_ncnst*this->nhl[0]];
this->iaddr2 = new int[2*this->nsp*this->nsp];
this->iaddr3 = new int[2*this->nsp*this->nsp*this->nsp];
for (int i = 0; i < 2*this->nsp*this->nsp; i++) {
this-> iaddr2[i] = -1;
for (int j = 0; j < this->nsp; j++) {
this->iaddr3[j+i*nsp] = -1;
}
}
//read in consts
int icmb[3] = {0}; // three is hardcoded maximum n in n-body potential
int iap = 0;
int jap = 0;
int kap = 0;
int nsf1 = 0, nsf2 = 0;
for (int i = 0; i < this->nhl[0]; i++) {
const_file >> this->itype[i];
for (int j = 0; i < this->ncomb_type[itype[i]]; i++) {
const_file >> icmb[j];
}
for (int j = 0; i < this->ncnst_type[itype[i]]; i++) {
const_file >> this->cnst[j+i*nhl[0]];
}
// distinguish n-body potentials
if (itype[i] < 100) {
if ((icmb[0] != iap) || (icmb[1] != jap)) {
this->iaddr2[0 + this->nsp*(icmb[0] + nsp*icmb[1])] = i+1;
this->iaddr2[0 + this->nsp*(icmb[1] + nsp*icmb[0])] = i+1;
}
this->iaddr2[1 + this->nsp*(icmb[0] + nsp*icmb[1])] = i+1;
this->iaddr2[1 + this->nsp*(icmb[1] + nsp*icmb[0])] = i+1;
nsf1++;
iap = icmb[0];
jap = icmb[1];
}
else if (itype[i] < 200){
if ((icmb[0] != iap) || (icmb[1] != jap) || (icmb[2] != kap)) {
this->iaddr3[0 + this->nsp*(icmb[0] + nsp*(icmb[1] + nsp*icmb[2]))] = i+1;
this->iaddr3[0 + this->nsp*(icmb[0] + nsp*(icmb[2] + nsp*icmb[1]))] = i+1;
}
this->iaddr3[1 + this->nsp*(icmb[0] + nsp*(icmb[1] + nsp*icmb[2]))] = i+1;
this->iaddr3[1 + this->nsp*(icmb[0] + nsp*(icmb[2] + nsp*icmb[1]))] = i+1;
nsf2++;
iap = icmb[0];
jap = icmb[1];
kap = icmb[2];
}
else {
error->all(FLERR,"Unknown itype");
}
}
if (this-> nhl[0] != nsf1 + nsf2) {
error->all(FLERR, "[Error] nsf.ne.nsf1+nsf2 !!!");
}
// read in weights
double nwgt[nl+1];
for (int i = 0 ; i < nl+1; i++) {
nwgt[i] = nhl[i] * nhl[i+1];
}
std::ifstream weights_file(arg[3]);
// TODO check file existence etc
int ncoeff;
weights_file >> ncoeff >> this->rcin >> this->rc3;
if (rc3 > rcin) {
rc3 = rcin;
// TODO add warning as in ryo_force_NN.F90:680
}
int nc = 0;
for (int i = 0; i < nl+1; i++) {
nc += nwgt[i];
}
if( ncoeff != nc ) {
error->all(FLERR, "[Error] num of parameters is not correct !!!");
}
if (nl == 1) {
this->wgt11 = new double[nhl[0] * nhl[1]];
this->wgt12 = new double[nhl[1]];
this->wgt21 = new double[1];
this->wgt22 = new double[1];
this->wgt23 = new double[1];
for (int ihl0 = 0; ihl0 < nhl[0]; ihl0++) {
for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
weights_file >> this->wgt11[ihl0 + ihl1*nhl[0]];
}
}
for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
weights_file >> this->wgt12[ihl1];
}
} else if (nl == 2) {
this->wgt11 = new double[1];
this->wgt12 = new double[1];
this->wgt21 = new double[nhl[0] * nhl[1]];
this->wgt22 = new double[nhl[1] * nhl[2]];
this->wgt23 = new double[nhl[2]];
for (int ihl0 = 0; ihl0 < nhl[0]; ihl0++) {
for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
weights_file >> this->wgt21[ihl0 + ihl1*nhl[0]];
}
}
for (int ihl0 = 0; ihl0 < nhl[1]; ihl0++) {
for (int ihl1 = 0; ihl1 < nhl[2]; ihl1++) {
weights_file >> this->wgt22[ihl0 + ihl1*nhl[0]];
}
}
for (int ihl1 = 0; ihl1 < nhl[2]; ihl1++) {
weights_file >> this->wgt23[ihl1];
}
}
// allocate arrays for weights
this->allocated = true;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairNeuroNet::init_style()
{
// Don't understand this yet, but seems to be what lammps wants here
int irequest = neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairNeuroNet::init_one(int i, int j)
{
return this->rcin;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairNeuroNet::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]) {
// fwrite(&epsilon[i][j],sizeof(double),1,fp);
// fwrite(&sigma[i][j],sizeof(double),1,fp);
// fwrite(&cut[i][j],sizeof(double),1,fp);
// }
// }
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairNeuroNet::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) {
// fread(&epsilon[i][j],sizeof(double),1,fp);
// fread(&sigma[i][j],sizeof(double),1,fp);
// fread(&cut[i][j],sizeof(double),1,fp);
// }
// MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
// MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
// MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
// }
// }
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairNeuroNet::write_restart_settings(FILE *fp)
{
// think this is supposed to be empty for us
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairNeuroNet::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
// think this is supposed to be empty for us
}
//MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
//MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
//MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
//MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
}
const int PairNeuroNet::nlmax = 2;
const int PairNeuroNet::max_ncnst = 2;
c4science · Help