Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/src/CNEURONET/pair_cneuronet.cpp
f9a4fb669e70lammm-devel

Last Change
Raw File

pair_cneuronet.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 "pair_cneuronet.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
#include "string.h"
#include "math_const.h"
#include <fstream>
#include <iostream>
#include <type_traits>
#include "CNEURONET/neigh_container.hh"
#include <cmath>
#include <stdexcept>
#include <numeric>
using namespace LAMMPS_NS;
using namespace MathConst;
using uint = PairCNeuroNet::uint;
using SymFunType = PairCNeuroNet::SymFunType;
using i_iterator = PairCNeuroNet::i_iterator;
using ij_iterator = PairCNeuroNet::ij_iterator;
using index_pair = PairCNeuroNet::index_pair;
using param_ptr = PairCNeuroNet::param_ptr;
using param_map = PairCNeuroNet::param_map;
using WeightType = PairCNeuroNet::WeightType;
using W_col_it = PairCNeuroNet::W_col_it;
using dfun_ret = PairCNeuroNet::dfun_ret;
static const int threeD = 3;
/* ---------------------------------------------------------------------- */
PairCNeuroNet::PairCNeuroNet(LAMMPS * lmp)
: Pair(lmp)
{}
/* ---------------------------------------------------------------------- */
PairCNeuroNet::~PairCNeuroNet()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
}
}
/* ---------------------------------------------------------------------- */
enum class PairCNeuroNet::SymFunType: int
{ Gaussian=1,
cosine=2,
polynomial=3,
Morse=4,
angular=101};
std::istream & operator>> ( std::istream & in, SymFunType & self) {
using en = std::underlying_type<SymFunType>::type;
in >> *reinterpret_cast<en*>(&self);
}
/* ---------------------------------------------------------------------- */
constexpr uint PairCNeuroNet::get_ncnst(const SymFunType & func_type) {
switch (func_type) {
case SymFunType::Gaussian:
return 2;
case SymFunType::cosine:
return 1;
case SymFunType::polynomial:
return 1;
case SymFunType::Morse:
return 2;
case SymFunType::angular:
return 1;
}
}
/* ---------------------------------------------------------------------- */
constexpr uint PairCNeuroNet::get_ncomb(const SymFunType & func_type) {
switch (func_type) {
case SymFunType::Gaussian:
return 2;
case SymFunType::cosine:
return 2;
case SymFunType::polynomial:
return 2;
case SymFunType::Morse:
return 2;
case SymFunType::angular:
return 3;
}
}
/* ---------------------------------------------------------------------- */
void PairCNeuroNet::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
size_t 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;
// Note: weights will be looped over along the number of
// input nodes, so that should be the contiguous dimension
this->weights.emplace_back(std::array<size_t, 2>({dummy, this->nhl.back()}));
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
this->weights.emplace_back(std::array<size_t, 2>({1, this->nhl.back()}));
nhl.push_back(1);
// 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->get_ncomb(itype[i]); j++) {
const_file >> icmb[j];
icmb[j]--;
}
// for now only for 2-body:
this->sym_fun_params[ordered(icmb[0]+1,icmb[1]+1)][itype[i]].push_back
({&this->cnst[i*this->max_ncnst], i});
double tmp;
for (int j = 0; j < this->get_ncnst(itype[i]); j++) {
const_file >> tmp;
this->cnst[j+i*this->max_ncnst] = tmp;
}
// distinguish n-body potentials
switch (this->get_ncnst(itype[i])) {
case 2 :
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];
break;
case 3 :
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];
break;
default:
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 !!!");
}
double fdummy1, fdummy2;
if (nl == 1) {
for (size_t ihl0 = 0; ihl0 < nhl[0]; ihl0++) {
for (size_t ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
weights_file >> this->weights[0].at({ihl1, ihl0}) >> fdummy1 >> fdummy2;
}
}
for (size_t ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
weights_file >> this->weights[1].at({0, ihl1}) >> fdummy1 >> fdummy2;
}
} else if (nl == 2) {
for (size_t ihl0 = 0; ihl0 < nhl[0]; ihl0++) {
for (size_t ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
weights_file >> this->weights[0][{ihl1, ihl0}] >> fdummy1 >> fdummy2;
}
}
for (size_t ihl0 = 0; ihl0 < nhl[1]; ihl0++) {
for (size_t ihl1 = 0; ihl1 < nhl[2]; ihl1++) {
weights_file >> this->weights[1][{ihl1, ihl0}] >> fdummy1 >> fdummy2;
}
}
for (size_t ihl1 = 0; ihl1 < nhl[2]; ihl1++) {
weights_file >> this->weights[2][{0, ihl1}] >> fdummy1 >> fdummy2;
}
}
// allocate arrays for weights
this->allocated = true;
}
/* ---------------------------------------------------------------------- */
inline double PairCNeuroNet::comp_r_ij(const double * pos_i, const double * pos_j) {
double r2{0.};
for (auto i = 0; i < threeD; ++i) {
auto delta = pos_i[i]-pos_j[i];
r2 += delta*delta;
}
return sqrt(r2);
}
/* ---------------------------------------------------------------------- */
inline double PairCNeuroNet::comp_fc_ij(const double & r_ij) {
if (r_ij < this->alpha * this->rcin) {
return 1.;
} else if (r_ij < this->rcin) {
return .5 * (
cos((r_ij-this->alpha*this->rcin)/
(this->rcin*(1.-this->alpha))*M_PI)
+ 1. );
} else {
return 0.;
}
}
/* ---------------------------------------------------------------------- */
inline double PairCNeuroNet::comp_Dfc_ij(const double & r_ij) {
if (r_ij < this->alpha * this->rcin) {
return 0.;
} else if (r_ij < this->rcin) {
return -.5 * M_PI/(this->rcin*(1.-this->alpha))*(
sin((r_ij-this->alpha*this->rcin)/
(this->rcin*(1.-this->alpha))*M_PI)
+ 1. );
} else {
return 0.;
}
}
/* ---------------------------------------------------------------------- */
template<SymFunType func_type>
inline double PairCNeuroNet::comp_fs_nij(const double * const param,
const double & r_ij) {
switch (func_type) {
case SymFunType::Gaussian: {
const auto & eta_n = param[0];
const auto delta_r = r_ij - param[1];
return exp(-eta_n*delta_r*delta_r);
}
default:
throw std::runtime_error("Function type not specified");
}
}
/* ---------------------------------------------------------------------- */
template<SymFunType func_type>
inline double PairCNeuroNet::comp_Dfs_nij(const double * const param,
const double & r_ij) {
switch (func_type) {
case SymFunType::Gaussian: {
const auto & eta_n = param[0];
const auto delta_r = r_ij - param[1];
return -2*eta_n*delta_r*exp(-eta_n*delta_r*delta_r);
}
default:
throw std::runtime_error("Function type not specified");
}
}
/* ---------------------------------------------------------------------- */
inline double PairCNeuroNet::comp_f_a(const double & x) {
return 1./(1.+exp(-x)) - .5;
}
/* ---------------------------------------------------------------------- */
inline double PairCNeuroNet::comp_Df_a(const double & x) {
const auto s = this->comp_f_a(x)+.5;
return s*(1-s);
}
/* ---------------------------------------------------------------------- */
template<SymFunType func_type>
inline void PairCNeuroNet::comp_G_imj(const ij_iterator & ij,
const param_ptr& param_vec) {
double fs_nij, Dfs_nij;
for (const auto & ptr: param_vec) {
const auto * param = ptr.first;
const auto & m2 = ptr.second;
fs_nij = this-> comp_fs_nij <func_type>(param, this->r_ij[ij]);
Dfs_nij = this->comp_Dfs_nij <func_type>(param, this->r_ij[ij]);
auto incr = fs_nij * this->fc_ij[ij];
this-> G_im (ij.get_j())[m2] += incr;
this-> G_im (ij.get_i())[m2] += incr;
this->DG_imj(ij )[m2] = Dfs_nij * this->fc_ij[ij] + fs_nij * this->Dfc_ij[ij];
}
}
/* ---------------------------------------------------------------------- */
inline double PairCNeuroNet::comp_ymm(W_row_it & w_it,
const i_iterator & at_i,
const size_t layer) {
auto TODO = comp_phi(w_it, at_i, layer);
return this->comp_f_a(TODO);
}
/* ---------------------------------------------------------------------- */
inline dfun_ret PairCNeuroNet::comp_Dymm(W_row_it & w_it,
const i_iterator & at_i,
const ij_iterator & at_k,
const size_t layer) {
auto && phi_Dphi = comp_Dphi(w_it, at_i, at_k, layer);
return {this->comp_f_a (phi_Dphi[0]),
this->comp_Df_a(phi_Dphi[0])*phi_Dphi[1]};
}
/* ---------------------------------------------------------------------- */
double PairCNeuroNet::comp_phi(W_row_it & w_it,
const i_iterator & at_i,
const size_t layer) {
switch (layer) {
case 0:
return std::inner_product(w_it.begin(), w_it.end(),
this->G_im(at_i), 0.);
break;
default:
auto && lower_layer = layer-1;
auto w_lower_it = this->weights[lower_layer].begin();
double phi{0.};
for (auto wmn_it = w_it.begin(); wmn_it != w_it.end(); ++wmn_it, ++w_lower_it) {
const auto &wmn =*wmn_it;
phi += wmn * this->comp_ymm(w_lower_it, at_i, layer-1);
}
return phi;
break;
}
}
/* ---------------------------------------------------------------------- */
dfun_ret PairCNeuroNet::comp_Dphi(W_row_it & w_it,
const i_iterator & at_i,
const ij_iterator & at_k,
const size_t layer) {
switch (layer) {
case 0: {
auto counter = 0;
double Gim{0}, DGim{0};
for (auto it_expl = w_it.begin(); it_expl != w_it.end(); ++it_expl, ++counter) {
Gim += (*it_expl)* G_im (at_i)[counter];
DGim += (*it_expl)*DG_imj(at_k)[counter];
}
// counter = 0;
// for (auto & w:w_it) {
// std::cout << counter++ << ", " << w << std::endl;
// }
auto wbeg = w_it.begin();
for (size_t i{0}; i < 120; ++i) {
auto wmn = *wbeg;
++wbeg;
}
dfun_ret bla = {std::inner_product(w_it.begin(), w_it.end(),
this->G_im (at_i), 0.),
std::inner_product(w_it.begin(), w_it.end(),
this->DG_imj(at_k), 0.)};
return bla;
break;
}
default:
auto && lower_layer = layer-1;
auto w_lower_it = this->weights[lower_layer].begin();
double phi{0.}, Dphi{0.};
for (auto wmn_it = w_it.begin(); wmn_it != w_it.end(); ++wmn_it, ++w_lower_it) {
const auto &wmn =*wmn_it;
auto && y_Dy = this->comp_Dymm(w_lower_it, at_i, at_k, lower_layer);
auto TODO = y_Dy;
phi += wmn * y_Dy[0];
Dphi += wmn * y_Dy[1];
}
return {phi, Dphi};
break;
}
}
/* ---------------------------------------------------------------------- */
template<int vflag_at, int eflag_at>
inline void PairCNeuroNet::comp_roll_um_sums_inner(){
for (auto & at_i: this->neighs) {
// innermost weights are last
auto w_row = this->weights.back().begin(); //iterator to the 1st and only
auto && E_iS = this->comp_phi(w_row, at_i, this->nl);
this->eng_vdwl += E_iS;
if (eflag_at) {
at_i.get_eatom() += E_iS;
}
for (auto & at_k: at_i) {
std::valarray<double> del_ik(threeD);
for (size_t i = 0; i < threeD; ++i) {
del_ik[i] = at_k.get_x()[i] - at_i.get_x()[i];
}
std::valarray<double> n_ik = del_ik/sqrt((del_ik*del_ik).sum());
/* ---------------------------------------------------------------------- */
auto && phi_Dphi = this->comp_Dphi(w_row, at_i, at_k, this->nl);
auto TODO = phi_Dphi;
auto f_ik = n_ik*phi_Dphi[1]; // force
for (size_t i = 0; i < threeD ; ++i) {
at_i.get_f()[i] += f_ik[i];
at_k.get_f()[i] -= f_ik[i];
if (vflag_at) {
auto * vat_i = at_i.get_vatom();
auto * vat_k = at_k.get_vatom();
const static size_t sigma[][threeD] = {{0, 3, 4},
{3, 1, 5},
{4, 5, 2}};
const static double factor[][threeD] = {{0.50, 0.25, 0.25},
{0.25, 0.50, 0.25},
{0.25, 0.25, 0.50}};
for (size_t j = 0; j < threeD; ++j) {
auto incr = factor [i][j] * f_ik[i] * del_ik[j];
vat_i[sigma[i][j]] += incr;
vat_k[sigma[i][j]] += incr;
}
}
}
}
}
}
/* ---------------------------------------------------------------------- */
template<int vflag_at>
inline void PairCNeuroNet::comp_roll_um_sums(){
if (this->eflag_atom) {
this->comp_roll_um_sums_inner<vflag_at, 1>();
} else {
this->comp_roll_um_sums_inner<vflag_at, 0>();
}
}
/* ---------------------------------------------------------------------- */
void PairCNeuroNet::compute(int eflag, int vflag) {
// lammps stuff
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
// the meat
this->neighs.set(this->list, this->atom, this->eatom, this->vatom);
this->r_ij.resize(this->neighs);
this-> fc_ij.resize(this->neighs);
this->Dfc_ij.resize(this->neighs);
// this->fs_nij.resize(this->neighs, this->nhl[0]);
// this->fs_Dnij.resize(this->neighs, this->nhl[0]);
this-> G_im .resize(this->neighs, this->nhl[0]);
this->DG_imj.resize(this->neighs, this->nhl[0]);
this-> G_im .set(0.);
this->DG_imj.set(0.);
// First, compute arrays that will be reused
for (auto & at_i:this->neighs) {
for (auto & ij:at_i) {
this->r_ij[ij] = comp_r_ij(at_i.get_x(), ij.get_x());
this->fc_ij[ij] = comp_fc_ij(this->r_ij[ij]);
this->Dfc_ij[ij] = comp_Dfc_ij(this->r_ij[ij]);
auto&& interaction = ordered(ij.get_type(), at_i.get_type());
size_t siz = sym_fun_params.size();
std::vector<std::array<int,2>> keys;
auto inter = interaction;
for(auto &item: sym_fun_params) keys.push_back(item.first);
for (auto & item: this->sym_fun_params.at(interaction)) {
const SymFunType & func_type = item.first;
const auto & param_vec = item.second;
switch (func_type) {
case SymFunType::Gaussian:
this->comp_G_imj<SymFunType::Gaussian>(ij, param_vec);
break;
default:
throw std::runtime_error("Function type not specified");
}
}
}
}
// now, roll up the sums
if (this->vflag_atom) {
this->comp_roll_um_sums<1>();
} else {
this->comp_roll_um_sums<0>();
}
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairCNeuroNet::settings(int narg, char **arg) // read_params
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairCNeuroNet::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=1;
neighbor->requests[irequest]->full=0;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairCNeuroNet::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);
}
using weigh_type = DimNumArray<2, double>;
template<>
void vector_write (const std::vector<weigh_type> & vec, FILE *fp) {
size_t siz = vec.size();
fwrite(&siz, sizeof(size_t), 1, fp);
for (const auto & weight: vec) {
weight.write(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);
fread(&vec[0], sizeof(T), siz, fp);
}
template<>
void vector_read(std::vector<weigh_type> & vec, FILE *fp) {
size_t siz;
fread(&siz, sizeof(size_t), 1, fp);
vec.resize(siz);
for(auto &weight:vec) {
weight.read(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 PairCNeuroNet::write_restart(FILE *fp)
{
write_restart_settings(fp);
scalar_write(this->nsp, fp);
scalar_write(this->nl, fp);
vector_write(this->weights, 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 PairCNeuroNet::read_restart(FILE *fp)
{
read_restart_settings(fp);
scalar_read(this->nsp, fp);
scalar_read(this->nl, fp);
vector_read(this->weights, 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 PairCNeuroNet::write_restart_settings(FILE *fp)
{
// think this is supposed to be empty for us
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCNeuroNet::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);
}
c4science · Help