Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86525840
pair_neuronet.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Oct 7, 00:54
Size
16 KB
Mime Type
text/x-c
Expires
Wed, Oct 9, 00:54 (2 d)
Engine
blob
Format
Raw Data
Handle
21439270
Attached To
rLAMMPS lammps
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
,
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
);
}
}
/* ---------------------------------------------------------------------- */
void
PairNeuroNet
::
compute
(
int
eflag
,
int
vflag
)
{
int
i
,
j
,
ii
,
jj
,
inum
,
jnum
,
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
->
nlocal
+
atom
->
nghost
;
//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
=
new
int
[(
nnmax
+
1
)
*
nmax
]
;
int
*
fort_type
=
new
int
[
nmax
]
;
double
*
fort_x
=
new
double
[
3
*
nmax
]
;
double
*
fort_f
=
new
double
[
3
*
nmax
]
;
double
*
fort_vatom
=
new
double
[
6
*
nmax
]
;
double
*
fort_eatom
=
new
double
[
nmax
]
;
double
*
hl1
=
new
double
[
nlocal
*
this
->
nhl
[
1
]]
;
double
*
hl2
=
new
double
[
nlocal
*
this
->
nhl
[
2
]]
;
double
*
gsf
=
new
double
[
nlocal
*
this
->
nhl
[
0
]]
;
double
*
dgsf
=
new
double
[
3
*
nhl
[
0
]
*
(
nmax
+
1
)
*
nlocal
]
;
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
]
+
1
;
}
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
;
int
fort_vflag_atom
=
bool
(
this
->
vflag_atom
);
//dummy = __nn_MOD_dsigmoid(&dummy);
//
__nn_MOD_force_nn
(
&
nmax
,
// namax
&
inum
,
// natm
fort_type
,
// tag
fort_x
,
// ra
&
nnmax
,
// nnmax
fort_f
,
// force
fort_vatom
,
// 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
,
// lspr
&
idummy
,
// mpi_world
&
idummy
,
// myid
fort_eatom
,
// 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
,
// hl1
hl2
,
// hl2
gsf
,
// gsf
dgsf
,
// 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
)
{
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
();
delete
[]
fort_neighbor
;
delete
[]
fort_type
;
delete
[]
fort_x
;
delete
[]
fort_f
;
delete
[]
fort_vatom
;
delete
[]
fort_eatom
;
delete
[]
hl1
;
delete
[]
hl2
;
delete
[]
gsf
;
delete
[]
dgsf
;
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
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
;
Event Timeline
Log In to Comment