Homec4science
Diffusion lammps (lammm-master-deprecated) Browse

lammps/src/NEURONET/pair_neuronet.cpp
a96d5dcf7eedlammm-master-deprecated

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>
#include <iostream>
#include<sys/types.h>
#include <unistd.h>
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, int * vflag_atom,
int * itype, double * cnst, int* iaddr2, int * iaddr3,
int * nsp, int * nhl, const int * max_ncnst, int * nl,
const int * nlmax, double * hl1, double *hl2,
double * gsf, double * dgsf, int * nal, int * nnl,
double * wgt11, double * wgt12, double * wgt21,
double * wgt22, double * wgt23, double * rc3);
}
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) {
memory->destroy(setflag);
memory->destroy(cutsq);
}
}
/* ---------------------------------------------------------------------- */
template<typename T>
inline void ensure_size(std::vector<T> & container, const size_t & size) {
if (container.capacity() < size) {
container.reserve(size);
}
}
void PairNeuroNet::compute(int eflag, int vflag)
{
double 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->nlocal + atom->nghost; //atom->nmax;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
int inum = list->inum; //number of atoms that have a neighbor list
int * ilist = list->ilist; //indices of those atoms
int * numneigh = list->numneigh; // number of their neighbors
int ** 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) {
int i = ilist[ii];
nnmax = (nnmax > numneigh[i]) ? nnmax : numneigh[i];
}
// ryo neighlist is matrix
ensure_size(this->fort_neighbor, (nnmax+1) * nmax);
for (int i = 0; i < nmax; i++) {
fort_neighbor[(nnmax+1) * i] = 0;
for (int j = 1; j < nnmax + 1; j++) {
fort_neighbor[(nnmax+1) * i + j] = -1;
}
}
ensure_size(this->fort_f, 3*nmax);
ensure_size(this->fort_vatom, 6*nmax);
ensure_size(this->fort_eatom,nmax);
ensure_size(this->hl1, nlocal * this->nhl[1]);
ensure_size(this->hl2, nlocal * this->nhl[2]);
ensure_size(this->gsf, nlocal * this->nhl[0]) ;
ensure_size(this->dgsf, static_cast<unsigned long int>(3)
* nhl[0] * static_cast<unsigned long int>(nnmax+1)
* nlocal) ;
for (int ii = 0; ii < inum; ++ii) {
int i = ilist[ii];
int *neigh_ptr = firstneigh[i];
int nb_neigh = numneigh[i];
fort_neighbor[0 + (nnmax+1)*i] = nb_neigh;
for (int j = 0; j < nb_neigh; j++) {
fort_neighbor[j+1 + (nnmax+1)*i] = (neigh_ptr[j] & NEIGHMASK) + 1;
}
}
int min = 0, max = 0;
for (int i = 0; i < (nnmax+1) * nmax; i++) {
min = (min < fort_neighbor[i]) ? min : fort_neighbor[i];
max = (max > fort_neighbor[i]) ? max : fort_neighbor[i];
}
// {
// int idumb = 0;
// char hostname[256];
// std::cout <<"PID " << ::getpid() << " ready for attach\n";
// fflush(stdout);
// while (0 == idumb)
// sleep(5);
// }
double h[] = {1,0,0,
0,1,0,
0,0,1};
double dummy = 0;
int idummy = 0;
double delta_eng_vdwl;
int fort_vflag_atom = bool(this->vflag_atom);
__nn_MOD_force_nn(&nmax, // namax
&inum, // natm
type, // tag
x[0], // ra
&nnmax, // nnmax
&fort_f[0], // force
&fort_vatom[0], // vatom
h, // h
h, // hi
&dummy, // tcom
&nghost, // nb
&idummy, // nbmax
&idummy, // lsb
&idummy, // nex
&idummy, // lsrc
&idummy, // myparity
&idummy, // nn
&dummy, // sv
&rcin, // rc
&fort_neighbor[0], // lspr
&idummy, // mpi_world
&idummy, // myid
&fort_eatom[0], // epi
&delta_eng_vdwl, // epot
&fort_vflag_atom, // vflag_atom
&itype[0], // itype
&cnst[0], // cnst
&iaddr2[0], // iaddr2
&iaddr3[0], // iaddr3
&nsp, // nsp
&nhl[0], // nhl
&max_ncnst, // max_ncnst
&nl, // nl
&nlmax, // nlmax
&hl1[0], // hl1
&hl2[0], // hl2
&gsf[0], // gsf
&dgsf[0], // dgsf
&nlocal, // nal
&nnmax, // nnl
&wgt11[0], // wgt11
&wgt12[0], // wgt12
&wgt21[0], // wgt21
&wgt22[0], // wgt22
&wgt23[0], // wgt23
&rc3); // rc3
this->eng_vdwl += delta_eng_vdwl;
// rewrite results into lammps form
for (int ii = 0; ii < inum; ++ii) {
int i = ilist[ii];
if (eflag_atom) {
eatom[i] += fort_eatom[ii];
}
for (int dir = 0; dir < 3; ++dir) {
f[i][dir] += fort_f[dir+3*ii];
}
if (vflag_atom) {
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 != 4)
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);
}
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
for (int i = 0; i < n+1; i++) {
for (int j = 0; j < n+1; j++) {
this->setflag[i][j] = 1;
}
}
// last layer always has one node
nhl.push_back(1);
while (nhl.size()<this->nlmax+2) {
nhl.push_back(0);
}
// allocate arrays for consts
this->itype.resize(this->nhl[0]);
this->cnst.resize(this->max_ncnst*this->nhl[0]);
this->iaddr2.resize(2*this->nsp*this->nsp);
this->iaddr3.resize(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 = -1;
int jap = -1;
int kap = -1;
int nsf1 = 0, nsf2 = 0;
for (int i = 0; i < this->nhl[0]; i++) {
const_file >> this->itype[i];
for (int j = 0; j < this->ncomb_type[itype[i]-1]; j++) {
const_file >> icmb[j];
icmb[j]--;
}
double tmp;
for (int j = 0; j < this->ncnst_type[itype[i]-1]; j++) {
const_file >> tmp;
this->cnst[j+i*this->max_ncnst] = tmp;
}
// distinguish n-body potentials
if (itype[i] < 100) {
if ((icmb[0] != iap) || (icmb[1] != jap)) {
this->iaddr2[0 + 2*(icmb[0] + nsp*icmb[1])] = i+1;
this->iaddr2[0 + 2*(icmb[1] + nsp*icmb[0])] = i+1;
}
this->iaddr2[1 + 2*(icmb[0] + nsp*icmb[1])] = i+1;
this->iaddr2[1 + 2*(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 + 2*(icmb[0] + nsp*(icmb[1] + nsp*icmb[2]))] = i+1;
this->iaddr3[0 + 2*(icmb[0] + nsp*(icmb[2] + nsp*icmb[1]))] = i+1;
}
this->iaddr3[1 + 2*(icmb[0] + nsp*(icmb[1] + nsp*icmb[2]))] = i+1;
this->iaddr3[1 + 2*(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
}
for (int i = 0; i < atom->ntypes+1; i++) {
for (int j = 0; j < atom->ntypes+1; j++) {
this->cutsq[i][j] = rcin*rcin;
}
}
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 !!!");
}
// allocate them anyways
this->wgt11.resize(nhl[0] * nhl[1]);
this->wgt12.resize(nhl[1]);
this->wgt21.resize(nhl[0] * nhl[1]);
this->wgt22.resize(nhl[1] * nhl[2]);
this->wgt23.resize(nhl[2]);
double fdummy1, fdummy2;
if (nl == 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]] >> fdummy1 >> fdummy2;
}
}
for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
weights_file >> this->wgt12[ihl1] >> fdummy1 >> fdummy2;
}
} else if (nl == 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]] >> fdummy1 >> fdummy2;
}
}
for (int ihl0 = 0; ihl0 < nhl[1]; ihl0++) {
for (int ihl1 = 0; ihl1 < nhl[2]; ihl1++) {
weights_file >> this->wgt22[ihl0 + ihl1*nhl[0]] >> fdummy1 >> fdummy2;
}
}
for (int ihl1 = 0; ihl1 < nhl[2]; ihl1++) {
weights_file >> this->wgt23[ihl1] >> fdummy1 >> fdummy2;
}
}
// 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);
neighbor->requests[irequest]->half=0;
neighbor->requests[irequest]->full=1;
}
/* ----------------------------------------------------------------------
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
------------------------------------------------------------------------- */
template<typename T>
void vector_write(const std::vector<T> & vec, FILE *fp) {
size_t siz = vec.size();
fwrite(&siz, sizeof(size_t), 1, fp);
fwrite(&vec[0], sizeof(T), siz, fp);
}
template<typename T>
void vector_read(std::vector<T> & vec, FILE *fp) {
size_t siz;
fread(&siz, sizeof(size_t), 1, fp);
vec.resize(siz);
fwrite(&vec[0], sizeof(T), siz, fp);
}
template<typename T>
void scalar_write(const T & scalar, FILE *fp) {
fwrite(&scalar, sizeof(T), 1, fp);
}
template<typename T>
void scalar_read(T & scalar, FILE *fp) {
fread(&scalar, sizeof(T), 1, fp);
}
void PairNeuroNet::write_restart(FILE *fp)
{
write_restart_settings(fp);
vector_write(this->wgt11, fp);
vector_write(this->wgt12, fp);
vector_write(this->wgt21, fp);
vector_write(this->wgt22, fp);
vector_write(this->wgt23, fp);
vector_write(this->nhl, fp);
vector_write(this->itype, fp);
vector_write(this->iaddr2, fp);
vector_write(this->iaddr3, fp);
vector_write(this->cnst, fp);
scalar_write(this->rcin, fp);
scalar_write(this->rc3, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairNeuroNet::read_restart(FILE *fp)
{
read_restart_settings(fp);
vector_read(this->wgt11, fp);
vector_read(this->wgt12, fp);
vector_read(this->wgt21, fp);
vector_read(this->wgt22, fp);
vector_read(this->wgt23, fp);
vector_read(this->nhl, fp);
vector_read(this->itype, fp);
vector_read(this->iaddr2, fp);
vector_read(this->iaddr3, fp);
vector_read(this->cnst, fp);
scalar_read(this->rcin, fp);
scalar_read(this->rc3, fp);
}
/* ----------------------------------------------------------------------
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