mirror of https://github.com/lammps/lammps.git
remove MEAM and REAX from lib folder
This commit is contained in:
Axel Kohlmeyer
parent
e6321e1020
commit
8b5887bfb1
|
|
@ -33,8 +33,6 @@ kokkos Kokkos package for GPU and many-core acceleration
|
|||
from Kokkos development team (Sandia)
|
||||
linalg set of BLAS and LAPACK routines needed by USER-ATC package
|
||||
from Axel Kohlmeyer (Temple U)
|
||||
meam modified embedded atom method (MEAM) potential, MEAM package
|
||||
from Greg Wagner (Sandia)
|
||||
message client/server communication library via MPI, sockets, files
|
||||
from Steve Plimpton (Sandia)
|
||||
molfile hooks to VMD molfile plugins, used by the USER-MOLFILE package
|
||||
|
|
@ -51,8 +49,6 @@ qmmm quantum mechanics/molecular mechanics coupling interface
|
|||
from Axel Kohlmeyer (Temple U)
|
||||
quip interface to QUIP/libAtoms framework, USER-QUIP package
|
||||
from Albert Bartok-Partay and Gabor Csanyi (U Cambridge)
|
||||
reax ReaxFF potential, REAX package
|
||||
from Adri van Duin (Penn State) and Aidan Thompson (Sandia)
|
||||
smd hooks to Eigen library, used by USER-SMD package
|
||||
from Georg Ganzenmueller (Ernst Mach Institute, Germany)
|
||||
voronoi hooks to the Voro++ library, used by compute voronoi/atom command
|
||||
|
|
|
|||
|
|
@ -1,9 +0,0 @@
|
|||
# dependencies. needed for parallel make
|
||||
$(DIR)meam_data.o: meam_data.F
|
||||
$(DIR)meam_cleanup.o: meam_cleanup.F $(DIR)meam_data.o
|
||||
$(DIR)meam_dens_final.o: meam_dens_final.F $(DIR)meam_data.o
|
||||
$(DIR)meam_dens_init.o: meam_dens_init.F $(DIR)meam_data.o
|
||||
$(DIR)meam_force.o: meam_force.F $(DIR)meam_data.o
|
||||
$(DIR)meam_setup_done.o: meam_setup_done.F $(DIR)meam_data.o
|
||||
$(DIR)meam_setup_global.o: meam_setup_global.F $(DIR)meam_data.o
|
||||
$(DIR)meam_setup_param.o: meam_setup_param.F $(DIR)meam_data.o
|
||||
|
|
@ -1 +0,0 @@
|
|||
*.mod
|
||||
|
|
@ -1 +0,0 @@
|
|||
../Install.py
|
||||
|
|
@ -1,57 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.gfortran
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F meam_cleanup.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o) fm_exp.o
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = g95
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
%.o:%.c
|
||||
$(CC) $(F90FLAGS) -c $<
|
||||
|
||||
include .depend
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
|
|
@ -1,61 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.gfortran
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F meam_cleanup.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o) fm_exp.o
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = gfortran
|
||||
CC = gcc
|
||||
F90FLAGS = -O3 -fPIC -ffast-math -ftree-vectorize -fexpensive-optimizations -fno-second-underscore -g
|
||||
#F90FLAGS = -O
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
LINK = g++
|
||||
LINKFLAGS = -O
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
%.o:%.c
|
||||
$(CC) $(F90FLAGS) -c $<
|
||||
|
||||
include .depend
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
|
|
@ -1,57 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.ifort
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F meam_cleanup.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o) fm_exp.o
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = ifort
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
%.o:%.c
|
||||
$(CC) $(F90FLAGS) -c $<
|
||||
|
||||
include .depend
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
# Settings that the LAMMPS build will import when this package library is used
|
||||
|
||||
meam_SYSINC =
|
||||
meam_SYSLIB =
|
||||
meam_SYSPATH =
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
# Settings that the LAMMPS build will import when this package library is used
|
||||
|
||||
meam_SYSINC =
|
||||
meam_SYSLIB = -lgfortran
|
||||
meam_SYSPATH =
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
# Settings that the LAMMPS build will import when this package library is used
|
||||
|
||||
meam_SYSINC =
|
||||
meam_SYSLIB = -lifcore -lsvml -lompstub -limf
|
||||
meam_SYSPATH = -L/opt/intel-11.1.046/lib/intel64
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
# Settings that the LAMMPS build will import when this package library is used
|
||||
|
||||
meam_SYSINC =
|
||||
meam_SYSLIB = -lifcore -lsvml -lompstub -limf
|
||||
meam_SYSPATH = -L/opt/intel/fce/10.0.023/lib
|
||||
|
|
@ -1,61 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.empty
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F meam_cleanup.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o) fm_exp.o
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = mpifort
|
||||
CC = mpicc
|
||||
F90FLAGS = -O3 -fPIC
|
||||
#F90FLAGS = -O
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
LINK = mpicxx
|
||||
LINKFLAGS = -O
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
%.o:%.c
|
||||
$(CC) $(F90FLAGS) -c $<
|
||||
|
||||
include .depend
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
|
|
@ -1,57 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.pgf90
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F meam_cleanup.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o) fm_exp.o
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = pgf90
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
%.o:%.c
|
||||
$(CC) $(F90FLAGS) -c $<
|
||||
|
||||
include .depend
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
|
|
@ -1 +0,0 @@
|
|||
Makefile.gfortran
|
||||
|
|
@ -1,59 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.glory
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F meam_cleanup.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o) fm_exp.o
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = mpif90
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
LINK = g++
|
||||
LINKFLAGS = -O
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
%.o:%.c
|
||||
$(CC) $(F90FLAGS) -c $<
|
||||
|
||||
include .depend
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
MEAM (modified embedded atom method) library
|
||||
|
||||
Greg Wagner, Sandia National Labs
|
||||
gjwagne at sandia.gov
|
||||
Jan 2007
|
||||
|
||||
This library is in implementation of the MEAM potential, specifically
|
||||
designed to work with LAMMPS.
|
||||
|
||||
-------------------------------------------------
|
||||
|
||||
This directory has source files to build a library that LAMMPS
|
||||
links against when using the MEAM package.
|
||||
|
||||
This library must be built with a F90 compiler, before LAMMPS is
|
||||
built, so LAMMPS can link against it.
|
||||
|
||||
You can type "make lib-meam" from the src directory to see help on how
|
||||
to build this library via make commands, or you can do the same thing
|
||||
by typing "python Install.py" from within this directory, or you can
|
||||
do it manually by following the instructions below.
|
||||
|
||||
Build the library using one of the provided Makefile.* files or create
|
||||
your own, specific to your compiler and system. For example:
|
||||
|
||||
make -f Makefile.gfortran
|
||||
|
||||
When you are done building this library, two files should
|
||||
exist in this directory:
|
||||
|
||||
libmeam.a the library LAMMPS will link against
|
||||
Makefile.lammps settings the LAMMPS Makefile will import
|
||||
|
||||
Makefile.lammps is created by the make command, by copying one of the
|
||||
Makefile.lammps.* files. See the EXTRAMAKE setting at the top of the
|
||||
Makefile.* files.
|
||||
|
||||
IMPORTANT: You must examine the final Makefile.lammps to insure it is
|
||||
correct for your system, else the LAMMPS build will likely fail.
|
||||
|
||||
Makefile.lammps has settings for 3 variables:
|
||||
|
||||
user-meam_SYSINC = leave blank for this package
|
||||
user-meam_SYSLIB = auxiliary F90 libs needed to link a F90 lib with
|
||||
a C++ program (LAMMPS) via a C++ compiler
|
||||
user-meam_SYSPATH = path(s) to where those libraries are
|
||||
|
||||
Because you have a F90 compiler on your system, you should have these
|
||||
libraries. But you will have to figure out which ones are needed and
|
||||
where they are. Examples of common configurations are in the
|
||||
Makefile.lammps.* files.
|
||||
|
|
@ -1,133 +0,0 @@
|
|||
/*
|
||||
Copyright (c) 2012,2013 Axel Kohlmeyer <akohlmey@gmail.com>
|
||||
All rights reserved.
|
||||
|
||||
Redistribution and use in source and binary forms, with or without
|
||||
modification, are permitted provided that the following conditions
|
||||
are met:
|
||||
|
||||
* Redistributions of source code must retain the above copyright
|
||||
notice, this list of conditions and the following disclaimer.
|
||||
* Redistributions in binary form must reproduce the above copyright
|
||||
notice, this list of conditions and the following disclaimer in the
|
||||
documentation and/or other materials provided with the distribution.
|
||||
* Neither the name of the <organization> nor the
|
||||
names of its contributors may be used to endorse or promote products
|
||||
derived from this software without specific prior written permission.
|
||||
|
||||
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
||||
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
ARE DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
|
||||
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
|
||||
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
|
||||
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
|
||||
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
||||
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
|
||||
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
||||
*/
|
||||
|
||||
/* faster versions of 2**x, e**x, and 10**x in single and double precision.
|
||||
*
|
||||
* Based on the Cephes math library 2.8
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
|
||||
/* internal definitions for the fastermath library */
|
||||
|
||||
/* IEEE 754 double precision floating point data manipulation */
|
||||
typedef union
|
||||
{
|
||||
double f;
|
||||
uint64_t u;
|
||||
struct {int32_t i0,i1;};
|
||||
} udi_t;
|
||||
#define FM_DOUBLE_BIAS 1023
|
||||
#define FM_DOUBLE_EMASK 2146435072
|
||||
#define FM_DOUBLE_MBITS 20
|
||||
#define FM_DOUBLE_MMASK 1048575
|
||||
#define FM_DOUBLE_EZERO 1072693248
|
||||
|
||||
/* generate 2**num in floating point by bitshifting */
|
||||
#define FM_DOUBLE_INIT_EXP(var,num) \
|
||||
var.i0 = 0; \
|
||||
var.i1 = (((int) num) + FM_DOUBLE_BIAS) << 20
|
||||
|
||||
/* double precision constants */
|
||||
#define FM_DOUBLE_LOG2OFE 1.4426950408889634074
|
||||
#define FM_DOUBLE_LOGEOF2 6.9314718055994530942e-1
|
||||
#define FM_DOUBLE_LOG2OF10 3.32192809488736234789
|
||||
#define FM_DOUBLE_LOG10OF2 3.0102999566398119521e-1
|
||||
#define FM_DOUBLE_LOG10OFE 4.3429448190325182765e-1
|
||||
#define FM_DOUBLE_SQRT2 1.41421356237309504880
|
||||
#define FM_DOUBLE_SQRTH 0.70710678118654752440
|
||||
|
||||
/* optimizer friendly implementation of exp2(x).
|
||||
*
|
||||
* strategy:
|
||||
*
|
||||
* split argument into an integer part and a fraction:
|
||||
* ipart = floor(x+0.5);
|
||||
* fpart = x - ipart;
|
||||
*
|
||||
* compute exp2(ipart) from setting the ieee754 exponent
|
||||
* compute exp2(fpart) using a pade' approximation for x in [-0.5;0.5[
|
||||
*
|
||||
* the result becomes: exp2(x) = exp2(ipart) * exp2(fpart)
|
||||
*/
|
||||
|
||||
static const double fm_exp2_q[] = {
|
||||
/* 1.00000000000000000000e0, */
|
||||
2.33184211722314911771e2,
|
||||
4.36821166879210612817e3
|
||||
};
|
||||
static const double fm_exp2_p[] = {
|
||||
2.30933477057345225087e-2,
|
||||
2.02020656693165307700e1,
|
||||
1.51390680115615096133e3
|
||||
};
|
||||
|
||||
static double fm_exp2(double x)
|
||||
{
|
||||
double ipart, fpart, px, qx;
|
||||
udi_t epart;
|
||||
|
||||
ipart = floor(x+0.5);
|
||||
fpart = x - ipart;
|
||||
FM_DOUBLE_INIT_EXP(epart,ipart);
|
||||
|
||||
x = fpart*fpart;
|
||||
|
||||
px = fm_exp2_p[0];
|
||||
px = px*x + fm_exp2_p[1];
|
||||
qx = x + fm_exp2_q[0];
|
||||
px = px*x + fm_exp2_p[2];
|
||||
qx = qx*x + fm_exp2_q[1];
|
||||
|
||||
px = px * fpart;
|
||||
|
||||
x = 1.0 + 2.0*(px/(qx-px));
|
||||
return epart.f*x;
|
||||
}
|
||||
|
||||
double fm_exp_(double *x)
|
||||
{
|
||||
#if defined(__BYTE_ORDER__)
|
||||
#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
|
||||
return fm_exp2(FM_DOUBLE_LOG2OFE * (*x));
|
||||
#endif
|
||||
#endif
|
||||
return exp(*x);
|
||||
}
|
||||
|
||||
/*
|
||||
* Local Variables:
|
||||
* mode: c
|
||||
* compile-command: "make -C .."
|
||||
* c-basic-offset: 4
|
||||
* fill-column: 76
|
||||
* indent-tabs-mode: nil
|
||||
* End:
|
||||
*/
|
||||
|
|
@ -1,26 +0,0 @@
|
|||
c Declaration in pair_meam.h:
|
||||
c
|
||||
c void meam_cleanup()
|
||||
c
|
||||
c Call from PairMEAM destructor
|
||||
c
|
||||
c meam_cleanup()
|
||||
c
|
||||
|
||||
subroutine meam_cleanup
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer dealloc_error
|
||||
|
||||
deallocate(phir,STAT=dealloc_error)
|
||||
deallocate(phirar,STAT=dealloc_error)
|
||||
deallocate(phirar1,STAT=dealloc_error)
|
||||
deallocate(phirar2,STAT=dealloc_error)
|
||||
deallocate(phirar3,STAT=dealloc_error)
|
||||
deallocate(phirar4,STAT=dealloc_error)
|
||||
deallocate(phirar5,STAT=dealloc_error)
|
||||
deallocate(phirar6,STAT=dealloc_error)
|
||||
|
||||
return
|
||||
end
|
||||
|
|
@ -1,87 +0,0 @@
|
|||
|
||||
module meam_data
|
||||
|
||||
integer, parameter :: maxelt = 5
|
||||
real*8 , external :: fm_exp
|
||||
|
||||
c cutforce = force cutoff
|
||||
c cutforcesq = force cutoff squared
|
||||
|
||||
real*8 cutforce,cutforcesq
|
||||
|
||||
c Ec_meam = cohesive energy
|
||||
c re_meam = nearest-neighbor distance
|
||||
c Omega_meam = atomic volume
|
||||
c B_meam = bulk modulus
|
||||
c Z_meam = number of first neighbors for reference structure
|
||||
c ielt_meam = atomic number of element
|
||||
c A_meam = adjustable parameter
|
||||
c alpha_meam = sqrt(9*Omega*B/Ec)
|
||||
c rho0_meam = density scaling parameter
|
||||
c delta_meam = heat of formation for alloys
|
||||
c beta[0-3]_meam = electron density constants
|
||||
c t[0-3]_meam = coefficients on densities in Gamma computation
|
||||
c rho_ref_meam = background density for reference structure
|
||||
c ibar_meam(i) = selection parameter for Gamma function for elt i,
|
||||
c lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j)
|
||||
c neltypes = maximum number of element type defined
|
||||
c eltind = index number of pair (similar to Voigt notation; ij = ji)
|
||||
c phir = pair potential function array
|
||||
c phirar[1-6] = spline coeffs
|
||||
c attrac_meam = attraction parameter in Rose energy
|
||||
c repuls_meam = repulsion parameter in Rose energy
|
||||
c nn2_meam = 1 if second nearest neighbors are to be computed, else 0
|
||||
c zbl_meam = 1 if zbl potential for small r to be use, else 0
|
||||
c emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0
|
||||
c bkgd_dyn = 1 if reference densities follows Dynamo, else 0
|
||||
c Cmin_meam, Cmax_meam = min and max values in screening cutoff
|
||||
c rc_meam = cutoff distance for meam
|
||||
c delr_meam = cutoff region for meam
|
||||
c ebound_meam = factor giving maximum boundary of sceen fcn ellipse
|
||||
c augt1 = flag for whether t1 coefficient should be augmented
|
||||
c ialloy = flag for newer alloy formulation (as in dynamo code)
|
||||
c mix_ref_t = flag to recover "old" way of computing t in reference config
|
||||
c erose_form = selection parameter for form of E_rose function
|
||||
c gsmooth_factor = factor determining length of G smoothing region
|
||||
c vind[23]D = Voight notation index maps for 2 and 3D
|
||||
c v2D,v3D = array of factors to apply for Voight notation
|
||||
|
||||
c nr,dr = pair function discretization parameters
|
||||
c nrar,rdrar = spline coeff array parameters
|
||||
|
||||
real*8 Ec_meam(maxelt,maxelt),re_meam(maxelt,maxelt)
|
||||
real*8 Omega_meam(maxelt),Z_meam(maxelt)
|
||||
real*8 A_meam(maxelt),alpha_meam(maxelt,maxelt),rho0_meam(maxelt)
|
||||
real*8 delta_meam(maxelt,maxelt)
|
||||
real*8 beta0_meam(maxelt),beta1_meam(maxelt)
|
||||
real*8 beta2_meam(maxelt),beta3_meam(maxelt)
|
||||
real*8 t0_meam(maxelt),t1_meam(maxelt)
|
||||
real*8 t2_meam(maxelt),t3_meam(maxelt)
|
||||
real*8 rho_ref_meam(maxelt)
|
||||
integer ibar_meam(maxelt),ielt_meam(maxelt)
|
||||
character*3 lattce_meam(maxelt,maxelt)
|
||||
integer nn2_meam(maxelt,maxelt)
|
||||
integer zbl_meam(maxelt,maxelt)
|
||||
integer eltind(maxelt,maxelt)
|
||||
integer neltypes
|
||||
|
||||
real*8, allocatable :: phir(:,:)
|
||||
|
||||
real*8, allocatable :: phirar(:,:),phirar1(:,:),phirar2(:,:),
|
||||
$ phirar3(:,:),phirar4(:,:),phirar5(:,:),phirar6(:,:)
|
||||
|
||||
real*8 attrac_meam(maxelt,maxelt),repuls_meam(maxelt,maxelt)
|
||||
|
||||
real*8 Cmin_meam(maxelt,maxelt,maxelt)
|
||||
real*8 Cmax_meam(maxelt,maxelt,maxelt)
|
||||
real*8 rc_meam,delr_meam,ebound_meam(maxelt,maxelt)
|
||||
integer augt1, ialloy, mix_ref_t, erose_form
|
||||
integer emb_lin_neg, bkgd_dyn
|
||||
real*8 gsmooth_factor
|
||||
integer vind2D(3,3),vind3D(3,3,3)
|
||||
integer v2D(6),v3D(10)
|
||||
|
||||
integer nr,nrar
|
||||
real*8 dr,rdrar
|
||||
|
||||
end module
|
||||
|
|
@ -1,296 +0,0 @@
|
|||
c Extern "C" declaration has the form:
|
||||
c
|
||||
c void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *,
|
||||
c int *, int *, int *,
|
||||
c double *, double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, double *, int *);
|
||||
c
|
||||
c Call from pair_meam.cpp has the form:
|
||||
c
|
||||
c meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
|
||||
c &eng_vdwl,eatom,ntype,type,fmap,
|
||||
c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
|
||||
c &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
|
||||
c dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
|
||||
c
|
||||
|
||||
subroutine meam_dens_final(nlocal, nmax,
|
||||
$ eflag_either, eflag_global, eflag_atom, eng_vdwl, eatom,
|
||||
$ ntype, type, fmap,
|
||||
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave,
|
||||
$ Gamma, dGamma1, dGamma2, dGamma3,
|
||||
$ rho, rho0, rho1, rho2, rho3, fp, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer nlocal, nmax, eflag_either, eflag_global, eflag_atom
|
||||
integer ntype, type, fmap
|
||||
real*8 eng_vdwl, eatom, Arho1, Arho2
|
||||
real*8 Arho2b, Arho3, Arho3b
|
||||
real*8 t_ave, tsq_ave
|
||||
real*8 Gamma, dGamma1, dGamma2, dGamma3
|
||||
real*8 rho, rho0, rho1, rho2, rho3
|
||||
real*8 fp
|
||||
integer errorflag
|
||||
|
||||
dimension eatom(nmax)
|
||||
dimension type(nmax), fmap(ntype)
|
||||
dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax)
|
||||
dimension Arho3(10,nmax), Arho3b(3,nmax), t_ave(3,nmax)
|
||||
dimension tsq_ave(3,nmax)
|
||||
dimension Gamma(nmax), dGamma1(nmax), dGamma2(nmax)
|
||||
dimension dGamma3(nmax), rho(nmax), rho0(nmax)
|
||||
dimension rho1(nmax), rho2(nmax), rho3(nmax)
|
||||
dimension fp(nmax)
|
||||
|
||||
integer i, elti
|
||||
integer m
|
||||
real*8 rhob, G, dG, Gbar, dGbar, gam, shp(3), shpi(3), Z
|
||||
real*8 B, denom, rho_bkgd
|
||||
|
||||
c Complete the calculation of density
|
||||
|
||||
do i = 1,nlocal
|
||||
|
||||
elti = fmap(type(i))
|
||||
if (elti.gt.0) then
|
||||
rho1(i) = 0.d0
|
||||
rho2(i) = -1.d0/3.d0*Arho2b(i)*Arho2b(i)
|
||||
rho3(i) = 0.d0
|
||||
do m = 1,3
|
||||
rho1(i) = rho1(i) + Arho1(m,i)*Arho1(m,i)
|
||||
rho3(i) = rho3(i) - 3.d0/5.d0*Arho3b(m,i)*Arho3b(m,i)
|
||||
enddo
|
||||
do m = 1,6
|
||||
rho2(i) = rho2(i) + v2D(m)*Arho2(m,i)*Arho2(m,i)
|
||||
enddo
|
||||
do m = 1,10
|
||||
rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i)
|
||||
enddo
|
||||
|
||||
if( rho0(i) .gt. 0.0 ) then
|
||||
if (ialloy.eq.1) then
|
||||
if (tsq_ave(1,i) .ne. 0.0d0) then
|
||||
t_ave(1,i) = t_ave(1,i)/tsq_ave(1,i)
|
||||
else
|
||||
t_ave(1,i) = 0.0d0
|
||||
endif
|
||||
if (tsq_ave(2,i) .ne. 0.0d0) then
|
||||
t_ave(2,i) = t_ave(2,i)/tsq_ave(2,i)
|
||||
else
|
||||
t_ave(2,i) = 0.0d0
|
||||
endif
|
||||
if (tsq_ave(3,i) .ne. 0.0d0) then
|
||||
t_ave(3,i) = t_ave(3,i)/tsq_ave(3,i)
|
||||
else
|
||||
t_ave(3,i) = 0.0d0
|
||||
endif
|
||||
else if (ialloy.eq.2) then
|
||||
t_ave(1,i) = t1_meam(elti)
|
||||
t_ave(2,i) = t2_meam(elti)
|
||||
t_ave(3,i) = t3_meam(elti)
|
||||
else
|
||||
t_ave(1,i) = t_ave(1,i)/rho0(i)
|
||||
t_ave(2,i) = t_ave(2,i)/rho0(i)
|
||||
t_ave(3,i) = t_ave(3,i)/rho0(i)
|
||||
endif
|
||||
endif
|
||||
|
||||
Gamma(i) = t_ave(1,i)*rho1(i)
|
||||
$ + t_ave(2,i)*rho2(i) + t_ave(3,i)*rho3(i)
|
||||
|
||||
if( rho0(i) .gt. 0.0 ) then
|
||||
Gamma(i) = Gamma(i)/(rho0(i)*rho0(i))
|
||||
end if
|
||||
|
||||
Z = Z_meam(elti)
|
||||
|
||||
call G_gam(Gamma(i),ibar_meam(elti),
|
||||
$ gsmooth_factor,G,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
call get_shpfcn(shp,lattce_meam(elti,elti))
|
||||
if (ibar_meam(elti).le.0) then
|
||||
Gbar = 1.d0
|
||||
dGbar = 0.d0
|
||||
else
|
||||
if (mix_ref_t.eq.1) then
|
||||
gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2)
|
||||
$ +t_ave(3,i)*shp(3))/(Z*Z)
|
||||
else
|
||||
gam = (t1_meam(elti)*shp(1)+t2_meam(elti)*shp(2)
|
||||
$ +t3_meam(elti)*shp(3))/(Z*Z)
|
||||
endif
|
||||
call G_gam(gam,ibar_meam(elti),gsmooth_factor,
|
||||
$ Gbar,errorflag)
|
||||
endif
|
||||
rho(i) = rho0(i) * G
|
||||
|
||||
if (mix_ref_t.eq.1) then
|
||||
if (ibar_meam(elti).le.0) then
|
||||
Gbar = 1.d0
|
||||
dGbar = 0.d0
|
||||
else
|
||||
gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2)
|
||||
$ +t_ave(3,i)*shp(3))/(Z*Z)
|
||||
call dG_gam(gam,ibar_meam(elti),gsmooth_factor,
|
||||
$ Gbar,dGbar)
|
||||
endif
|
||||
rho_bkgd = rho0_meam(elti)*Z*Gbar
|
||||
else
|
||||
if (bkgd_dyn.eq.1) then
|
||||
rho_bkgd = rho0_meam(elti)*Z
|
||||
else
|
||||
rho_bkgd = rho_ref_meam(elti)
|
||||
endif
|
||||
endif
|
||||
rhob = rho(i)/rho_bkgd
|
||||
denom = 1.d0/rho_bkgd
|
||||
|
||||
call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG)
|
||||
|
||||
dGamma1(i) = (G - 2*dG*Gamma(i))*denom
|
||||
|
||||
if( rho0(i) .ne. 0.d0 ) then
|
||||
dGamma2(i) = (dG/rho0(i))*denom
|
||||
else
|
||||
dGamma2(i) = 0.d0
|
||||
end if
|
||||
|
||||
c dGamma3 is nonzero only if we are using the "mixed" rule for
|
||||
c computing t in the reference system (which is not correct, but
|
||||
c included for backward compatibility
|
||||
if (mix_ref_t.eq.1) then
|
||||
dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom
|
||||
else
|
||||
dGamma3(i) = 0.0
|
||||
endif
|
||||
|
||||
B = A_meam(elti)*Ec_meam(elti,elti)
|
||||
|
||||
if( rhob .ne. 0.d0 ) then
|
||||
if (emb_lin_neg.eq.1 .and. rhob.le.0) then
|
||||
fp(i) = -B
|
||||
else
|
||||
fp(i) = B*(log(rhob)+1.d0)
|
||||
endif
|
||||
if (eflag_either.ne.0) then
|
||||
if (eflag_global.ne.0) then
|
||||
if (emb_lin_neg.eq.1 .and. rhob.le.0) then
|
||||
eng_vdwl = eng_vdwl - B*rhob
|
||||
else
|
||||
eng_vdwl = eng_vdwl + B*rhob*log(rhob)
|
||||
endif
|
||||
endif
|
||||
if (eflag_atom.ne.0) then
|
||||
if (emb_lin_neg.eq.1 .and. rhob.le.0) then
|
||||
eatom(i) = eatom(i) - B*rhob
|
||||
else
|
||||
eatom(i) = eatom(i) + B*rhob*log(rhob)
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
else
|
||||
if (emb_lin_neg.eq.1) then
|
||||
fp(i) = -B
|
||||
else
|
||||
fp(i) = B
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine G_gam(Gamma,ibar,gsmooth_factor,G,errorflag)
|
||||
c Compute G(Gamma) based on selection flag ibar:
|
||||
c 0 => G = sqrt(1+Gamma)
|
||||
c 1 => G = exp(Gamma/2)
|
||||
c 2 => not implemented
|
||||
c 3 => G = 2/(1+exp(-Gamma))
|
||||
c 4 => G = sqrt(1+Gamma)
|
||||
c -5 => G = +-sqrt(abs(1+Gamma))
|
||||
use meam_data , only: fm_exp
|
||||
implicit none
|
||||
real*8 Gamma,G
|
||||
real*8 gsmooth_factor, gsmooth_switchpoint
|
||||
integer ibar, errorflag
|
||||
if (ibar.eq.0.or.ibar.eq.4) then
|
||||
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1)
|
||||
if (Gamma.lt.gsmooth_switchpoint) then
|
||||
c e.g. gsmooth_factor is 99, then:
|
||||
c gsmooth_switchpoint = -0.99
|
||||
c G = 0.01*(-0.99/Gamma)**99
|
||||
G = 1/(gsmooth_factor+1)
|
||||
$ *(gsmooth_switchpoint/Gamma)**gsmooth_factor
|
||||
G = sqrt(G)
|
||||
else
|
||||
G = sqrt(1.d0+Gamma)
|
||||
endif
|
||||
else if (ibar.eq.1) then
|
||||
G = fm_exp(Gamma/2.d0)
|
||||
else if (ibar.eq.3) then
|
||||
G = 2.d0/(1.d0+exp(-Gamma))
|
||||
else if (ibar.eq.-5) then
|
||||
if ((1.d0+Gamma).ge.0) then
|
||||
G = sqrt(1.d0+Gamma)
|
||||
else
|
||||
G = -sqrt(-1.d0-Gamma)
|
||||
endif
|
||||
else
|
||||
errorflag = 1
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dG_gam(Gamma,ibar,gsmooth_factor,G,dG)
|
||||
c Compute G(Gamma) and dG(gamma) based on selection flag ibar:
|
||||
c 0 => G = sqrt(1+Gamma)
|
||||
c 1 => G = fm_exp(Gamma/2)
|
||||
c 2 => not implemented
|
||||
c 3 => G = 2/(1+fm_exp(-Gamma))
|
||||
c 4 => G = sqrt(1+Gamma)
|
||||
c -5 => G = +-sqrt(abs(1+Gamma))
|
||||
use meam_data , only: fm_exp
|
||||
real*8 Gamma,G,dG
|
||||
real*8 gsmooth_factor, gsmooth_switchpoint
|
||||
integer ibar
|
||||
if (ibar.eq.0.or.ibar.eq.4) then
|
||||
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1)
|
||||
if (Gamma.lt.gsmooth_switchpoint) then
|
||||
c e.g. gsmooth_factor is 99, then:
|
||||
c gsmooth_switchpoint = -0.99
|
||||
c G = 0.01*(-0.99/Gamma)**99
|
||||
G = 1/(gsmooth_factor+1)
|
||||
$ *(gsmooth_switchpoint/Gamma)**gsmooth_factor
|
||||
G = sqrt(G)
|
||||
dG = -gsmooth_factor*G/(2.0*Gamma)
|
||||
else
|
||||
G = sqrt(1.d0+Gamma)
|
||||
dG = 1.d0/(2.d0*G)
|
||||
endif
|
||||
else if (ibar.eq.1) then
|
||||
G = fm_exp(Gamma/2.d0)
|
||||
dG = G/2.d0
|
||||
else if (ibar.eq.3) then
|
||||
G = 2.d0/(1.d0+fm_exp(-Gamma))
|
||||
dG = G*(2.d0-G)/2
|
||||
else if (ibar.eq.-5) then
|
||||
if ((1.d0+Gamma).ge.0) then
|
||||
G = sqrt(1.d0+Gamma)
|
||||
dG = 1.d0/(2.d0*G)
|
||||
else
|
||||
G = -sqrt(-1.d0-Gamma)
|
||||
dG = -1.d0/(2.d0*G)
|
||||
endif
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
|
@ -1,564 +0,0 @@
|
|||
c Extern "C" declaration has the form:
|
||||
c
|
||||
c void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *,
|
||||
c int *, int *, int *, int *,
|
||||
c double *, double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, double *, int *);
|
||||
c
|
||||
c
|
||||
c Call from pair_meam.cpp has the form:
|
||||
c
|
||||
c meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0],
|
||||
c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],
|
||||
c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
|
||||
c rho0,&arho1[0][0],&arho2[0][0],arho2b,
|
||||
c &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag);
|
||||
c
|
||||
|
||||
subroutine meam_dens_init(i, nmax,
|
||||
$ ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ numneigh_full, firstneigh_full,
|
||||
$ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b,
|
||||
$ arho3, arho3b, t_ave, tsq_ave, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer i, nmax, ntype, type, fmap
|
||||
real*8 x
|
||||
integer numneigh, firstneigh, numneigh_full, firstneigh_full
|
||||
real*8 scrfcn, dscrfcn, fcpair
|
||||
real*8 rho0, arho1, arho2
|
||||
real*8 arho2b, arho3, arho3b, t_ave, tsq_ave
|
||||
integer errorflag
|
||||
integer j,jn
|
||||
|
||||
dimension x(3,nmax)
|
||||
dimension type(nmax), fmap(ntype)
|
||||
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
|
||||
dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh)
|
||||
dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax)
|
||||
dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax)
|
||||
dimension t_ave(3,nmax), tsq_ave(3,nmax)
|
||||
|
||||
errorflag = 0
|
||||
|
||||
c Compute screening function and derivatives
|
||||
call getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ numneigh_full, firstneigh_full,
|
||||
$ ntype, type, fmap)
|
||||
|
||||
c Calculate intermediate density terms to be communicated
|
||||
call calc_rho1(i, nmax, ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
|
||||
$ arho3, arho3b, t_ave, tsq_ave)
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ numneigh_full, firstneigh_full,
|
||||
$ ntype, type, fmap)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer i, nmax
|
||||
real*8 scrfcn, dscrfcn, fcpair, x
|
||||
integer numneigh, firstneigh, numneigh_full, firstneigh_full
|
||||
integer ntype, type, fmap
|
||||
|
||||
dimension scrfcn(numneigh), dscrfcn(numneigh)
|
||||
dimension fcpair(numneigh), x(3,nmax)
|
||||
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
|
||||
dimension type(nmax), fmap(ntype)
|
||||
|
||||
integer jn,j,kn,k
|
||||
integer elti,eltj,eltk
|
||||
real*8 xitmp,yitmp,zitmp,delxij,delyij,delzij,rij2,rij
|
||||
real*8 xjtmp,yjtmp,zjtmp,delxik,delyik,delzik,rik2,rik
|
||||
real*8 xktmp,yktmp,zktmp,delxjk,delyjk,delzjk,rjk2,rjk
|
||||
real*8 xik,xjk,sij,fcij,sfcij,dfcij,sikj,dfikj,cikj
|
||||
real*8 Cmin,Cmax,delc,ebound,rbound,a,coef1,coef2
|
||||
real*8 coef1a,coef1b,coef2a,coef2b
|
||||
real*8 dcikj
|
||||
real*8 dC1a,dC1b,dC2a,dC2b
|
||||
real*8 rnorm,fc,dfc,drinv
|
||||
|
||||
drinv = 1.d0/delr_meam
|
||||
elti = fmap(type(i))
|
||||
|
||||
if (elti.gt.0) then
|
||||
|
||||
xitmp = x(1,i)
|
||||
yitmp = x(2,i)
|
||||
zitmp = x(3,i)
|
||||
|
||||
do jn = 1,numneigh
|
||||
j = firstneigh(jn)
|
||||
|
||||
eltj = fmap(type(j))
|
||||
if (eltj.gt.0) then
|
||||
|
||||
c First compute screening function itself, sij
|
||||
xjtmp = x(1,j)
|
||||
yjtmp = x(2,j)
|
||||
zjtmp = x(3,j)
|
||||
delxij = xjtmp - xitmp
|
||||
delyij = yjtmp - yitmp
|
||||
delzij = zjtmp - zitmp
|
||||
rij2 = delxij*delxij + delyij*delyij + delzij*delzij
|
||||
rij = sqrt(rij2)
|
||||
if (rij.gt.rc_meam) then
|
||||
fcij = 0.0
|
||||
dfcij = 0.d0
|
||||
sij = 0.d0
|
||||
else
|
||||
rnorm = (rc_meam-rij)*drinv
|
||||
call screen(i, j, nmax, x, rij2, sij,
|
||||
$ numneigh_full, firstneigh_full, ntype, type, fmap)
|
||||
call dfcut(rnorm,fc,dfc)
|
||||
fcij = fc
|
||||
dfcij = dfc*drinv
|
||||
endif
|
||||
|
||||
c Now compute derivatives
|
||||
dscrfcn(jn) = 0.d0
|
||||
sfcij = sij*fcij
|
||||
if (sfcij.eq.0.d0.or.sfcij.eq.1.d0) goto 100
|
||||
rbound = ebound_meam(elti,eltj) * rij2
|
||||
do kn = 1,numneigh_full
|
||||
k = firstneigh_full(kn)
|
||||
if (k.eq.j) goto 10
|
||||
eltk = fmap(type(k))
|
||||
if (eltk.eq.0) goto 10
|
||||
xktmp = x(1,k)
|
||||
yktmp = x(2,k)
|
||||
zktmp = x(3,k)
|
||||
delxjk = xktmp - xjtmp
|
||||
delyjk = yktmp - yjtmp
|
||||
delzjk = zktmp - zjtmp
|
||||
rjk2 = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk
|
||||
if (rjk2.gt.rbound) goto 10
|
||||
delxik = xktmp - xitmp
|
||||
delyik = yktmp - yitmp
|
||||
delzik = zktmp - zitmp
|
||||
rik2 = delxik*delxik + delyik*delyik + delzik*delzik
|
||||
if (rik2.gt.rbound) goto 10
|
||||
xik = rik2/rij2
|
||||
xjk = rjk2/rij2
|
||||
a = 1 - (xik-xjk)*(xik-xjk)
|
||||
c if a < 0, then ellipse equation doesn't describe this case and
|
||||
c atom k can't possibly screen i-j
|
||||
if (a.le.0.d0) goto 10
|
||||
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
|
||||
Cmax = Cmax_meam(elti,eltj,eltk)
|
||||
Cmin = Cmin_meam(elti,eltj,eltk)
|
||||
if (cikj.ge.Cmax) then
|
||||
goto 10
|
||||
c Note that cikj may be slightly negative (within numerical
|
||||
c tolerance) if atoms are colinear, so don't reject that case here
|
||||
c (other negative cikj cases were handled by the test on "a" above)
|
||||
c Note that we never have 0<cikj<Cmin here, else sij=0 (rejected above)
|
||||
else
|
||||
delc = Cmax - Cmin
|
||||
cikj = (cikj-Cmin)/delc
|
||||
call dfcut(cikj,sikj,dfikj)
|
||||
coef1 = dfikj/(delc*sikj)
|
||||
call dCfunc(rij2,rik2,rjk2,dCikj)
|
||||
dscrfcn(jn) = dscrfcn(jn) + coef1*dCikj
|
||||
endif
|
||||
10 continue
|
||||
enddo
|
||||
coef1 = sfcij
|
||||
coef2 = sij*dfcij/rij
|
||||
dscrfcn(jn) = dscrfcn(jn)*coef1 - coef2
|
||||
100 continue
|
||||
|
||||
scrfcn(jn) = sij
|
||||
fcpair(jn) = fcij
|
||||
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine calc_rho1(i, nmax, ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
|
||||
$ arho3, arho3b, t_ave, tsq_ave)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer i, nmax, ntype, type, fmap
|
||||
real*8 x
|
||||
integer numneigh, firstneigh
|
||||
real*8 scrfcn, fcpair, rho0, arho1, arho2
|
||||
real*8 arho2b, arho3, arho3b, t_ave, tsq_ave
|
||||
|
||||
dimension type(nmax), fmap(ntype), x(3,nmax)
|
||||
dimension firstneigh(numneigh)
|
||||
dimension scrfcn(numneigh), fcpair(numneigh)
|
||||
dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax)
|
||||
dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax)
|
||||
dimension t_ave(3,nmax), tsq_ave(3,nmax)
|
||||
|
||||
integer jn,j,m,n,p,elti,eltj
|
||||
integer nv2,nv3
|
||||
real*8 xtmp,ytmp,ztmp,delij(3),rij2,rij,sij
|
||||
real*8 ai,aj,rhoa0j,rhoa1j,rhoa2j,rhoa3j,A1j,A2j,A3j
|
||||
real*8 G,Gbar,gam,shp(3)
|
||||
real*8 ro0i,ro0j
|
||||
real*8 rhoa0i,rhoa1i,rhoa2i,rhoa3i,A1i,A2i,A3i
|
||||
|
||||
elti = fmap(type(i))
|
||||
xtmp = x(1,i)
|
||||
ytmp = x(2,i)
|
||||
ztmp = x(3,i)
|
||||
do jn = 1,numneigh
|
||||
if (scrfcn(jn).ne.0.d0) then
|
||||
j = firstneigh(jn)
|
||||
sij = scrfcn(jn)*fcpair(jn)
|
||||
delij(1) = x(1,j) - xtmp
|
||||
delij(2) = x(2,j) - ytmp
|
||||
delij(3) = x(3,j) - ztmp
|
||||
rij2 = delij(1)*delij(1) + delij(2)*delij(2)
|
||||
$ + delij(3)*delij(3)
|
||||
if (rij2.lt.cutforcesq) then
|
||||
eltj = fmap(type(j))
|
||||
rij = sqrt(rij2)
|
||||
ai = rij/re_meam(elti,elti) - 1.d0
|
||||
aj = rij/re_meam(eltj,eltj) - 1.d0
|
||||
ro0i = rho0_meam(elti)
|
||||
ro0j = rho0_meam(eltj)
|
||||
rhoa0j = ro0j*fm_exp(-beta0_meam(eltj)*aj)*sij
|
||||
rhoa1j = ro0j*fm_exp(-beta1_meam(eltj)*aj)*sij
|
||||
rhoa2j = ro0j*fm_exp(-beta2_meam(eltj)*aj)*sij
|
||||
rhoa3j = ro0j*fm_exp(-beta3_meam(eltj)*aj)*sij
|
||||
rhoa0i = ro0i*fm_exp(-beta0_meam(elti)*ai)*sij
|
||||
rhoa1i = ro0i*fm_exp(-beta1_meam(elti)*ai)*sij
|
||||
rhoa2i = ro0i*fm_exp(-beta2_meam(elti)*ai)*sij
|
||||
rhoa3i = ro0i*fm_exp(-beta3_meam(elti)*ai)*sij
|
||||
if (ialloy.eq.1) then
|
||||
rhoa1j = rhoa1j * t1_meam(eltj)
|
||||
rhoa2j = rhoa2j * t2_meam(eltj)
|
||||
rhoa3j = rhoa3j * t3_meam(eltj)
|
||||
rhoa1i = rhoa1i * t1_meam(elti)
|
||||
rhoa2i = rhoa2i * t2_meam(elti)
|
||||
rhoa3i = rhoa3i * t3_meam(elti)
|
||||
endif
|
||||
rho0(i) = rho0(i) + rhoa0j
|
||||
rho0(j) = rho0(j) + rhoa0i
|
||||
c For ialloy = 2, use single-element value (not average)
|
||||
if (ialloy.ne.2) then
|
||||
t_ave(1,i) = t_ave(1,i) + t1_meam(eltj)*rhoa0j
|
||||
t_ave(2,i) = t_ave(2,i) + t2_meam(eltj)*rhoa0j
|
||||
t_ave(3,i) = t_ave(3,i) + t3_meam(eltj)*rhoa0j
|
||||
t_ave(1,j) = t_ave(1,j) + t1_meam(elti)*rhoa0i
|
||||
t_ave(2,j) = t_ave(2,j) + t2_meam(elti)*rhoa0i
|
||||
t_ave(3,j) = t_ave(3,j) + t3_meam(elti)*rhoa0i
|
||||
endif
|
||||
if (ialloy.eq.1) then
|
||||
tsq_ave(1,i) = tsq_ave(1,i) +
|
||||
$ t1_meam(eltj)*t1_meam(eltj)*rhoa0j
|
||||
tsq_ave(2,i) = tsq_ave(2,i) +
|
||||
$ t2_meam(eltj)*t2_meam(eltj)*rhoa0j
|
||||
tsq_ave(3,i) = tsq_ave(3,i) +
|
||||
$ t3_meam(eltj)*t3_meam(eltj)*rhoa0j
|
||||
tsq_ave(1,j) = tsq_ave(1,j) +
|
||||
$ t1_meam(elti)*t1_meam(elti)*rhoa0i
|
||||
tsq_ave(2,j) = tsq_ave(2,j) +
|
||||
$ t2_meam(elti)*t2_meam(elti)*rhoa0i
|
||||
tsq_ave(3,j) = tsq_ave(3,j) +
|
||||
$ t3_meam(elti)*t3_meam(elti)*rhoa0i
|
||||
endif
|
||||
Arho2b(i) = Arho2b(i) + rhoa2j
|
||||
Arho2b(j) = Arho2b(j) + rhoa2i
|
||||
|
||||
A1j = rhoa1j/rij
|
||||
A2j = rhoa2j/rij2
|
||||
A3j = rhoa3j/(rij2*rij)
|
||||
A1i = rhoa1i/rij
|
||||
A2i = rhoa2i/rij2
|
||||
A3i = rhoa3i/(rij2*rij)
|
||||
nv2 = 1
|
||||
nv3 = 1
|
||||
do m = 1,3
|
||||
Arho1(m,i) = Arho1(m,i) + A1j*delij(m)
|
||||
Arho1(m,j) = Arho1(m,j) - A1i*delij(m)
|
||||
Arho3b(m,i) = Arho3b(m,i) + rhoa3j*delij(m)/rij
|
||||
Arho3b(m,j) = Arho3b(m,j) - rhoa3i*delij(m)/rij
|
||||
do n = m,3
|
||||
Arho2(nv2,i) = Arho2(nv2,i) + A2j*delij(m)*delij(n)
|
||||
Arho2(nv2,j) = Arho2(nv2,j) + A2i*delij(m)*delij(n)
|
||||
nv2 = nv2+1
|
||||
do p = n,3
|
||||
Arho3(nv3,i) = Arho3(nv3,i)
|
||||
$ + A3j*delij(m)*delij(n)*delij(p)
|
||||
Arho3(nv3,j) = Arho3(nv3,j)
|
||||
$ - A3i*delij(m)*delij(n)*delij(p)
|
||||
nv3 = nv3+1
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine screen(i, j, nmax, x, rijsq, sij,
|
||||
$ numneigh_full, firstneigh_full, ntype, type, fmap)
|
||||
c Screening function
|
||||
c Inputs: i = atom 1 id (integer)
|
||||
c j = atom 2 id (integer)
|
||||
c rijsq = squared distance between i and j
|
||||
c Outputs: sij = screening function
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer i,j,nmax,k,nk,m
|
||||
real*8 x,rijsq,sij
|
||||
integer numneigh_full, firstneigh_full
|
||||
integer ntype, type, fmap
|
||||
|
||||
dimension x(3,nmax), firstneigh_full(numneigh_full)
|
||||
dimension type(nmax), fmap(ntype)
|
||||
|
||||
integer elti,eltj,eltk
|
||||
real*8 delxik,delyik,delzik
|
||||
real*8 delxjk,delyjk,delzjk
|
||||
real*8 riksq,rjksq,xik,xjk,cikj,a,delc,sikj,fcij,rij
|
||||
real*8 Cmax,Cmin,rbound
|
||||
|
||||
sij = 1.d0
|
||||
elti = fmap(type(i))
|
||||
eltj = fmap(type(j))
|
||||
|
||||
c if rjksq > ebound*rijsq, atom k is definitely outside the ellipse
|
||||
rbound = ebound_meam(elti,eltj)*rijsq
|
||||
|
||||
do nk = 1,numneigh_full
|
||||
k = firstneigh_full(nk)
|
||||
eltk = fmap(type(k))
|
||||
if (k.eq.j) goto 10
|
||||
delxjk = x(1,k) - x(1,j)
|
||||
delyjk = x(2,k) - x(2,j)
|
||||
delzjk = x(3,k) - x(3,j)
|
||||
rjksq = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk
|
||||
if (rjksq.gt.rbound) goto 10
|
||||
delxik = x(1,k) - x(1,i)
|
||||
delyik = x(2,k) - x(2,i)
|
||||
delzik = x(3,k) - x(3,i)
|
||||
riksq = delxik*delxik + delyik*delyik + delzik*delzik
|
||||
if (riksq.gt.rbound) goto 10
|
||||
xik = riksq/rijsq
|
||||
xjk = rjksq/rijsq
|
||||
a = 1 - (xik-xjk)*(xik-xjk)
|
||||
c if a < 0, then ellipse equation doesn't describe this case and
|
||||
c atom k can't possibly screen i-j
|
||||
if (a.le.0.d0) goto 10
|
||||
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
|
||||
Cmax = Cmax_meam(elti,eltj,eltk)
|
||||
Cmin = Cmin_meam(elti,eltj,eltk)
|
||||
if (cikj.ge.Cmax) then
|
||||
goto 10
|
||||
c note that cikj may be slightly negative (within numerical
|
||||
c tolerance) if atoms are colinear, so don't reject that case here
|
||||
c (other negative cikj cases were handled by the test on "a" above)
|
||||
else if (cikj.le.Cmin) then
|
||||
sij = 0.d0
|
||||
goto 20
|
||||
else
|
||||
delc = Cmax - Cmin
|
||||
cikj = (cikj-Cmin)/delc
|
||||
call fcut(cikj,sikj)
|
||||
endif
|
||||
sij = sij * sikj
|
||||
10 continue
|
||||
enddo
|
||||
|
||||
20 continue
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2,
|
||||
$ ntype,type,fmap,x,scrfcn,fcpair)
|
||||
c Inputs: i,j,k = id's of 3 atom triplet
|
||||
c jn = id of i-j pair
|
||||
c rij2 = squared distance between i and j
|
||||
c Outputs: dsij1 = deriv. of sij w.r.t. rik
|
||||
c dsij2 = deriv. of sij w.r.t. rjk
|
||||
use meam_data
|
||||
implicit none
|
||||
integer i,j,k,jn,nmax,numneigh
|
||||
integer elti,eltj,eltk
|
||||
real*8 rij2,rik2,rjk2,dsij1,dsij2
|
||||
integer ntype, type, fmap
|
||||
real*8 x, scrfcn, fcpair
|
||||
|
||||
dimension type(nmax), fmap(ntype)
|
||||
dimension x(3,nmax), scrfcn(numneigh), fcpair(numneigh)
|
||||
|
||||
real*8 dxik,dyik,dzik
|
||||
real*8 dxjk,dyjk,dzjk
|
||||
real*8 rbound,delc,sij,xik,xjk,cikj,sikj,dfc,a
|
||||
real*8 Cmax,Cmin,dCikj1,dCikj2
|
||||
|
||||
sij = scrfcn(jn)*fcpair(jn)
|
||||
elti = fmap(type(i))
|
||||
eltj = fmap(type(j))
|
||||
eltk = fmap(type(k))
|
||||
Cmax = Cmax_meam(elti,eltj,eltk)
|
||||
Cmin = Cmin_meam(elti,eltj,eltk)
|
||||
|
||||
dsij1 = 0.d0
|
||||
dsij2 = 0.d0
|
||||
if ((sij.ne.0.d0).and.(sij.ne.1.d0)) then
|
||||
rbound = rij2*ebound_meam(elti,eltj)
|
||||
delc = Cmax-Cmin
|
||||
dxjk = x(1,k) - x(1,j)
|
||||
dyjk = x(2,k) - x(2,j)
|
||||
dzjk = x(3,k) - x(3,j)
|
||||
rjk2 = dxjk*dxjk + dyjk*dyjk + dzjk*dzjk
|
||||
if (rjk2.le.rbound) then
|
||||
dxik = x(1,k) - x(1,i)
|
||||
dyik = x(2,k) - x(2,i)
|
||||
dzik = x(3,k) - x(3,i)
|
||||
rik2 = dxik*dxik + dyik*dyik + dzik*dzik
|
||||
if (rik2.le.rbound) then
|
||||
xik = rik2/rij2
|
||||
xjk = rjk2/rij2
|
||||
a = 1 - (xik-xjk)*(xik-xjk)
|
||||
if (a.ne.0.d0) then
|
||||
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
|
||||
if (cikj.ge.Cmin.and.cikj.le.Cmax) then
|
||||
cikj = (cikj-Cmin)/delc
|
||||
call dfcut(cikj,sikj,dfc)
|
||||
call dCfunc2(rij2,rik2,rjk2,dCikj1,dCikj2)
|
||||
a = sij/delc*dfc/sikj
|
||||
dsij1 = a*dCikj1
|
||||
dsij2 = a*dCikj2
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine fcut(xi,fc)
|
||||
c cutoff function
|
||||
implicit none
|
||||
real*8 xi,fc
|
||||
real*8 a
|
||||
if (xi.ge.1.d0) then
|
||||
fc = 1.d0
|
||||
else if (xi.le.0.d0) then
|
||||
fc = 0.d0
|
||||
else
|
||||
a = 1.d0-xi
|
||||
a = a*a
|
||||
a = a*a
|
||||
a = 1.d0-a
|
||||
fc = a*a
|
||||
c fc = xi
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dfcut(xi,fc,dfc)
|
||||
c cutoff function and its derivative
|
||||
implicit none
|
||||
real*8 xi,fc,dfc,a,a3,a4
|
||||
if (xi.ge.1.d0) then
|
||||
fc = 1.d0
|
||||
dfc = 0.d0
|
||||
else if (xi.le.0.d0) then
|
||||
fc = 0.d0
|
||||
dfc = 0.d0
|
||||
else
|
||||
a = 1.d0-xi
|
||||
a3 = a*a*a
|
||||
a4 = a*a3
|
||||
fc = (1.d0-a4)**2
|
||||
dfc = 8*(1.d0-a4)*a3
|
||||
c fc = xi
|
||||
c dfc = 1.d0
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dCfunc(rij2,rik2,rjk2,dCikj)
|
||||
c Inputs: rij,rij2,rik2,rjk2
|
||||
c Outputs: dCikj = derivative of Cikj w.r.t. rij
|
||||
implicit none
|
||||
real*8 rij2,rik2,rjk2,dCikj
|
||||
real*8 rij4,a,b,denom
|
||||
|
||||
rij4 = rij2*rij2
|
||||
a = rik2-rjk2
|
||||
b = rik2+rjk2
|
||||
denom = rij4 - a*a
|
||||
denom = denom*denom
|
||||
dCikj = -4*(-2*rij2*a*a + rij4*b + a*a*b)/denom
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dCfunc2(rij2,rik2,rjk2,dCikj1,dCikj2)
|
||||
c Inputs: rij,rij2,rik2,rjk2
|
||||
c Outputs: dCikj1 = derivative of Cikj w.r.t. rik
|
||||
c dCikj2 = derivative of Cikj w.r.t. rjk
|
||||
implicit none
|
||||
real*8 rij2,rik2,rjk2,dCikj1,dCikj2
|
||||
real*8 rij4,rik4,rjk4,a,b,denom
|
||||
|
||||
rij4 = rij2*rij2
|
||||
rik4 = rik2*rik2
|
||||
rjk4 = rjk2*rjk2
|
||||
a = rik2-rjk2
|
||||
b = rik2+rjk2
|
||||
denom = rij4 - a*a
|
||||
denom = denom*denom
|
||||
dCikj1 = 4*rij2*(rij4 + rik4 + 2*rik2*rjk2 - 3*rjk4 - 2*rij2*a)/
|
||||
$ denom
|
||||
dCikj2 = 4*rij2*(rij4 - 3*rik4 + 2*rik2*rjk2 + rjk4 + 2*rij2*a)/
|
||||
$ denom
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -1,608 +0,0 @@
|
|||
c Extern "C" declaration has the form:
|
||||
c
|
||||
c void meam_force_(int *, int *, int *, double *, int *, int *, int *, double *,
|
||||
c int *, int *, int *, int *, double *, double *,
|
||||
c double *, double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, double *, double *, int *);
|
||||
c
|
||||
c Call from pair_meam.cpp has the form:
|
||||
c
|
||||
c meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom,
|
||||
c &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
|
||||
c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],
|
||||
c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
|
||||
c dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
|
||||
c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
|
||||
c &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag);
|
||||
c
|
||||
|
||||
subroutine meam_force(i, nmax,
|
||||
$ eflag_either, eflag_global, eflag_atom, vflag_atom,
|
||||
$ eng_vdwl, eatom, ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh, numneigh_full, firstneigh_full,
|
||||
$ scrfcn, dscrfcn, fcpair,
|
||||
$ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp,
|
||||
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, f,
|
||||
$ vatom, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer eflag_either, eflag_global, eflag_atom, vflag_atom
|
||||
integer nmax, ntype, type, fmap
|
||||
real*8 eng_vdwl, eatom, x
|
||||
integer numneigh, firstneigh, numneigh_full, firstneigh_full
|
||||
real*8 scrfcn, dscrfcn, fcpair
|
||||
real*8 dGamma1, dGamma2, dGamma3
|
||||
real*8 rho0, rho1, rho2, rho3, fp
|
||||
real*8 Arho1, Arho2, Arho2b
|
||||
real*8 Arho3, Arho3b
|
||||
real*8 t_ave, tsq_ave, f, vatom
|
||||
integer errorflag
|
||||
|
||||
dimension eatom(nmax)
|
||||
dimension type(nmax), fmap(ntype)
|
||||
dimension x(3,nmax)
|
||||
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
|
||||
dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh)
|
||||
dimension dGamma1(nmax), dGamma2(nmax), dGamma3(nmax)
|
||||
dimension rho0(nmax), rho1(nmax), rho2(nmax), rho3(nmax), fp(nmax)
|
||||
dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax)
|
||||
dimension Arho3(10,nmax), Arho3b(3,nmax)
|
||||
dimension t_ave(3,nmax), tsq_ave(3,nmax), f(3,nmax), vatom(6,nmax)
|
||||
|
||||
integer i,j,jn,k,kn,kk,m,n,p,q
|
||||
integer nv2,nv3,elti,eltj,eltk,ind
|
||||
real*8 xitmp,yitmp,zitmp,delij(3),delref(3),rij2,rij,rij3
|
||||
real*8 delik(3),deljk(3),v(6),fi(3),fj(3)
|
||||
real*8 Eu,astar,astarp,third,sixth
|
||||
real*8 pp,phiforce,dUdrij,dUdsij,dUdrijm(3),force,forcem
|
||||
real*8 B,r,recip,phi,phip,rhop,a
|
||||
real*8 sij,fcij,dfcij,ds(3)
|
||||
real*8 a0,a1,a1i,a1j,a2,a2i,a2j
|
||||
real*8 a3i,a3j,a3i1,a3i2,a3j1,a3j2
|
||||
real*8 G,dG,Gbar,dGbar,gam,shpi(3),shpj(3),Z,denom
|
||||
real*8 ai,aj,ro0i,ro0j,invrei,invrej
|
||||
real*8 b0,rhoa0j,drhoa0j,rhoa0i,drhoa0i
|
||||
real*8 b1,rhoa1j,drhoa1j,rhoa1i,drhoa1i
|
||||
real*8 b2,rhoa2j,drhoa2j,rhoa2i,drhoa2i
|
||||
real*8 a3,a3a,b3,rhoa3j,drhoa3j,rhoa3i,drhoa3i
|
||||
real*8 drho0dr1,drho0dr2,drho0ds1,drho0ds2
|
||||
real*8 drho1dr1,drho1dr2,drho1ds1,drho1ds2
|
||||
real*8 drho1drm1(3),drho1drm2(3)
|
||||
real*8 drho2dr1,drho2dr2,drho2ds1,drho2ds2
|
||||
real*8 drho2drm1(3),drho2drm2(3)
|
||||
real*8 drho3dr1,drho3dr2,drho3ds1,drho3ds2
|
||||
real*8 drho3drm1(3),drho3drm2(3)
|
||||
real*8 dt1dr1,dt1dr2,dt1ds1,dt1ds2
|
||||
real*8 dt2dr1,dt2dr2,dt2ds1,dt2ds2
|
||||
real*8 dt3dr1,dt3dr2,dt3ds1,dt3ds2
|
||||
real*8 drhodr1,drhodr2,drhods1,drhods2,drhodrm1(3),drhodrm2(3)
|
||||
real*8 arg,arg1,arg2
|
||||
real*8 arg1i1,arg1j1,arg1i2,arg1j2,arg2i2,arg2j2
|
||||
real*8 arg1i3,arg1j3,arg2i3,arg2j3,arg3i3,arg3j3
|
||||
real*8 dsij1,dsij2,force1,force2
|
||||
real*8 t1i,t2i,t3i,t1j,t2j,t3j
|
||||
|
||||
errorflag = 0
|
||||
third = 1.0/3.0
|
||||
sixth = 1.0/6.0
|
||||
|
||||
c Compute forces atom i
|
||||
|
||||
elti = fmap(type(i))
|
||||
|
||||
if (elti.gt.0) then
|
||||
xitmp = x(1,i)
|
||||
yitmp = x(2,i)
|
||||
zitmp = x(3,i)
|
||||
|
||||
c Treat each pair
|
||||
do jn = 1,numneigh
|
||||
|
||||
j = firstneigh(jn)
|
||||
eltj = fmap(type(j))
|
||||
|
||||
if (scrfcn(jn).ne.0.d0.and.eltj.gt.0) then
|
||||
|
||||
sij = scrfcn(jn)*fcpair(jn)
|
||||
delij(1) = x(1,j) - xitmp
|
||||
delij(2) = x(2,j) - yitmp
|
||||
delij(3) = x(3,j) - zitmp
|
||||
rij2 = delij(1)*delij(1) + delij(2)*delij(2)
|
||||
$ + delij(3)*delij(3)
|
||||
if (rij2.lt.cutforcesq) then
|
||||
rij = sqrt(rij2)
|
||||
r = rij
|
||||
|
||||
c Compute phi and phip
|
||||
ind = eltind(elti,eltj)
|
||||
pp = rij*rdrar + 1.0D0
|
||||
kk = pp
|
||||
kk = min(kk,nrar-1)
|
||||
pp = pp - kk
|
||||
pp = min(pp,1.0D0)
|
||||
phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp
|
||||
$ + phirar1(kk,ind))*pp + phirar(kk,ind)
|
||||
phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp
|
||||
$ + phirar4(kk,ind)
|
||||
recip = 1.0d0/r
|
||||
|
||||
if (eflag_either.ne.0) then
|
||||
if (eflag_global.ne.0) eng_vdwl = eng_vdwl + phi*sij
|
||||
if (eflag_atom.ne.0) then
|
||||
eatom(i) = eatom(i) + 0.5*phi*sij
|
||||
eatom(j) = eatom(j) + 0.5*phi*sij
|
||||
endif
|
||||
endif
|
||||
|
||||
c write(1,*) "force_meamf: phi: ",phi
|
||||
c write(1,*) "force_meamf: phip: ",phip
|
||||
|
||||
c Compute pair densities and derivatives
|
||||
invrei = 1.d0/re_meam(elti,elti)
|
||||
ai = rij*invrei - 1.d0
|
||||
ro0i = rho0_meam(elti)
|
||||
rhoa0i = ro0i*fm_exp(-beta0_meam(elti)*ai)
|
||||
drhoa0i = -beta0_meam(elti)*invrei*rhoa0i
|
||||
rhoa1i = ro0i*fm_exp(-beta1_meam(elti)*ai)
|
||||
drhoa1i = -beta1_meam(elti)*invrei*rhoa1i
|
||||
rhoa2i = ro0i*fm_exp(-beta2_meam(elti)*ai)
|
||||
drhoa2i = -beta2_meam(elti)*invrei*rhoa2i
|
||||
rhoa3i = ro0i*fm_exp(-beta3_meam(elti)*ai)
|
||||
drhoa3i = -beta3_meam(elti)*invrei*rhoa3i
|
||||
|
||||
if (elti.ne.eltj) then
|
||||
invrej = 1.d0/re_meam(eltj,eltj)
|
||||
aj = rij*invrej - 1.d0
|
||||
ro0j = rho0_meam(eltj)
|
||||
rhoa0j = ro0j*fm_exp(-beta0_meam(eltj)*aj)
|
||||
drhoa0j = -beta0_meam(eltj)*invrej*rhoa0j
|
||||
rhoa1j = ro0j*fm_exp(-beta1_meam(eltj)*aj)
|
||||
drhoa1j = -beta1_meam(eltj)*invrej*rhoa1j
|
||||
rhoa2j = ro0j*fm_exp(-beta2_meam(eltj)*aj)
|
||||
drhoa2j = -beta2_meam(eltj)*invrej*rhoa2j
|
||||
rhoa3j = ro0j*fm_exp(-beta3_meam(eltj)*aj)
|
||||
drhoa3j = -beta3_meam(eltj)*invrej*rhoa3j
|
||||
else
|
||||
rhoa0j = rhoa0i
|
||||
drhoa0j = drhoa0i
|
||||
rhoa1j = rhoa1i
|
||||
drhoa1j = drhoa1i
|
||||
rhoa2j = rhoa2i
|
||||
drhoa2j = drhoa2i
|
||||
rhoa3j = rhoa3i
|
||||
drhoa3j = drhoa3i
|
||||
endif
|
||||
|
||||
if (ialloy.eq.1) then
|
||||
rhoa1j = rhoa1j * t1_meam(eltj)
|
||||
rhoa2j = rhoa2j * t2_meam(eltj)
|
||||
rhoa3j = rhoa3j * t3_meam(eltj)
|
||||
rhoa1i = rhoa1i * t1_meam(elti)
|
||||
rhoa2i = rhoa2i * t2_meam(elti)
|
||||
rhoa3i = rhoa3i * t3_meam(elti)
|
||||
drhoa1j = drhoa1j * t1_meam(eltj)
|
||||
drhoa2j = drhoa2j * t2_meam(eltj)
|
||||
drhoa3j = drhoa3j * t3_meam(eltj)
|
||||
drhoa1i = drhoa1i * t1_meam(elti)
|
||||
drhoa2i = drhoa2i * t2_meam(elti)
|
||||
drhoa3i = drhoa3i * t3_meam(elti)
|
||||
endif
|
||||
|
||||
nv2 = 1
|
||||
nv3 = 1
|
||||
arg1i1 = 0.d0
|
||||
arg1j1 = 0.d0
|
||||
arg1i2 = 0.d0
|
||||
arg1j2 = 0.d0
|
||||
arg1i3 = 0.d0
|
||||
arg1j3 = 0.d0
|
||||
arg3i3 = 0.d0
|
||||
arg3j3 = 0.d0
|
||||
do n = 1,3
|
||||
do p = n,3
|
||||
do q = p,3
|
||||
arg = delij(n)*delij(p)*delij(q)*v3D(nv3)
|
||||
arg1i3 = arg1i3 + Arho3(nv3,i)*arg
|
||||
arg1j3 = arg1j3 - Arho3(nv3,j)*arg
|
||||
nv3 = nv3+1
|
||||
enddo
|
||||
arg = delij(n)*delij(p)*v2D(nv2)
|
||||
arg1i2 = arg1i2 + Arho2(nv2,i)*arg
|
||||
arg1j2 = arg1j2 + Arho2(nv2,j)*arg
|
||||
nv2 = nv2+1
|
||||
enddo
|
||||
arg1i1 = arg1i1 + Arho1(n,i)*delij(n)
|
||||
arg1j1 = arg1j1 - Arho1(n,j)*delij(n)
|
||||
arg3i3 = arg3i3 + Arho3b(n,i)*delij(n)
|
||||
arg3j3 = arg3j3 - Arho3b(n,j)*delij(n)
|
||||
enddo
|
||||
|
||||
c rho0 terms
|
||||
drho0dr1 = drhoa0j * sij
|
||||
drho0dr2 = drhoa0i * sij
|
||||
|
||||
c rho1 terms
|
||||
a1 = 2*sij/rij
|
||||
drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1
|
||||
drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1
|
||||
a1 = 2.d0*sij/rij
|
||||
do m = 1,3
|
||||
drho1drm1(m) = a1*rhoa1j*Arho1(m,i)
|
||||
drho1drm2(m) = -a1*rhoa1i*Arho1(m,j)
|
||||
enddo
|
||||
|
||||
c rho2 terms
|
||||
a2 = 2*sij/rij2
|
||||
drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2
|
||||
$ - 2.d0/3.d0*Arho2b(i)*drhoa2j*sij
|
||||
drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2
|
||||
$ - 2.d0/3.d0*Arho2b(j)*drhoa2i*sij
|
||||
a2 = 4*sij/rij2
|
||||
do m = 1,3
|
||||
drho2drm1(m) = 0.d0
|
||||
drho2drm2(m) = 0.d0
|
||||
do n = 1,3
|
||||
drho2drm1(m) = drho2drm1(m)
|
||||
$ + Arho2(vind2D(m,n),i)*delij(n)
|
||||
drho2drm2(m) = drho2drm2(m)
|
||||
$ - Arho2(vind2D(m,n),j)*delij(n)
|
||||
enddo
|
||||
drho2drm1(m) = a2*rhoa2j*drho2drm1(m)
|
||||
drho2drm2(m) = -a2*rhoa2i*drho2drm2(m)
|
||||
enddo
|
||||
|
||||
c rho3 terms
|
||||
rij3 = rij*rij2
|
||||
a3 = 2*sij/rij3
|
||||
a3a = 6.d0/5.d0*sij/rij
|
||||
drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3
|
||||
$ - a3a*(drhoa3j - rhoa3j/rij)*arg3i3
|
||||
drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3
|
||||
$ - a3a*(drhoa3i - rhoa3i/rij)*arg3j3
|
||||
a3 = 6*sij/rij3
|
||||
a3a = 6*sij/(5*rij)
|
||||
do m = 1,3
|
||||
drho3drm1(m) = 0.d0
|
||||
drho3drm2(m) = 0.d0
|
||||
nv2 = 1
|
||||
do n = 1,3
|
||||
do p = n,3
|
||||
arg = delij(n)*delij(p)*v2D(nv2)
|
||||
drho3drm1(m) = drho3drm1(m)
|
||||
$ + Arho3(vind3D(m,n,p),i)*arg
|
||||
drho3drm2(m) = drho3drm2(m)
|
||||
$ + Arho3(vind3D(m,n,p),j)*arg
|
||||
nv2 = nv2 + 1
|
||||
enddo
|
||||
enddo
|
||||
drho3drm1(m) = (a3*drho3drm1(m) - a3a*Arho3b(m,i))
|
||||
$ *rhoa3j
|
||||
drho3drm2(m) = (-a3*drho3drm2(m) + a3a*Arho3b(m,j))
|
||||
$ *rhoa3i
|
||||
enddo
|
||||
|
||||
c Compute derivatives of weighting functions t wrt rij
|
||||
t1i = t_ave(1,i)
|
||||
t2i = t_ave(2,i)
|
||||
t3i = t_ave(3,i)
|
||||
t1j = t_ave(1,j)
|
||||
t2j = t_ave(2,j)
|
||||
t3j = t_ave(3,j)
|
||||
|
||||
if (ialloy.eq.1) then
|
||||
|
||||
a1i = 0.d0
|
||||
a1j = 0.d0
|
||||
a2i = 0.d0
|
||||
a2j = 0.d0
|
||||
a3i = 0.d0
|
||||
a3j = 0.d0
|
||||
if ( tsq_ave(1,i) .ne. 0.d0 ) then
|
||||
a1i = drhoa0j*sij/tsq_ave(1,i)
|
||||
endif
|
||||
if ( tsq_ave(1,j) .ne. 0.d0 ) then
|
||||
a1j = drhoa0i*sij/tsq_ave(1,j)
|
||||
endif
|
||||
if ( tsq_ave(2,i) .ne. 0.d0 ) then
|
||||
a2i = drhoa0j*sij/tsq_ave(2,i)
|
||||
endif
|
||||
if ( tsq_ave(2,j) .ne. 0.d0 ) then
|
||||
a2j = drhoa0i*sij/tsq_ave(2,j)
|
||||
endif
|
||||
if ( tsq_ave(3,i) .ne. 0.d0 ) then
|
||||
a3i = drhoa0j*sij/tsq_ave(3,i)
|
||||
endif
|
||||
if ( tsq_ave(3,j) .ne. 0.d0 ) then
|
||||
a3j = drhoa0i*sij/tsq_ave(3,j)
|
||||
endif
|
||||
|
||||
dt1dr1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2)
|
||||
dt1dr2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2)
|
||||
dt2dr1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2)
|
||||
dt2dr2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2)
|
||||
dt3dr1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2)
|
||||
dt3dr2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2)
|
||||
|
||||
else if (ialloy.eq.2) then
|
||||
|
||||
dt1dr1 = 0.d0
|
||||
dt1dr2 = 0.d0
|
||||
dt2dr1 = 0.d0
|
||||
dt2dr2 = 0.d0
|
||||
dt3dr1 = 0.d0
|
||||
dt3dr2 = 0.d0
|
||||
|
||||
else
|
||||
|
||||
ai = 0.d0
|
||||
if( rho0(i) .ne. 0.d0 ) then
|
||||
ai = drhoa0j*sij/rho0(i)
|
||||
end if
|
||||
aj = 0.d0
|
||||
if( rho0(j) .ne. 0.d0 ) then
|
||||
aj = drhoa0i*sij/rho0(j)
|
||||
end if
|
||||
|
||||
dt1dr1 = ai*(t1_meam(eltj)-t1i)
|
||||
dt1dr2 = aj*(t1_meam(elti)-t1j)
|
||||
dt2dr1 = ai*(t2_meam(eltj)-t2i)
|
||||
dt2dr2 = aj*(t2_meam(elti)-t2j)
|
||||
dt3dr1 = ai*(t3_meam(eltj)-t3i)
|
||||
dt3dr2 = aj*(t3_meam(elti)-t3j)
|
||||
|
||||
endif
|
||||
|
||||
c Compute derivatives of total density wrt rij, sij and rij(3)
|
||||
call get_shpfcn(shpi,lattce_meam(elti,elti))
|
||||
call get_shpfcn(shpj,lattce_meam(eltj,eltj))
|
||||
drhodr1 = dGamma1(i)*drho0dr1
|
||||
$ + dGamma2(i)*
|
||||
$ (dt1dr1*rho1(i)+t1i*drho1dr1
|
||||
$ + dt2dr1*rho2(i)+t2i*drho2dr1
|
||||
$ + dt3dr1*rho3(i)+t3i*drho3dr1)
|
||||
$ - dGamma3(i)*
|
||||
$ (shpi(1)*dt1dr1+shpi(2)*dt2dr1+shpi(3)*dt3dr1)
|
||||
drhodr2 = dGamma1(j)*drho0dr2
|
||||
$ + dGamma2(j)*
|
||||
$ (dt1dr2*rho1(j)+t1j*drho1dr2
|
||||
$ + dt2dr2*rho2(j)+t2j*drho2dr2
|
||||
$ + dt3dr2*rho3(j)+t3j*drho3dr2)
|
||||
$ - dGamma3(j)*
|
||||
$ (shpj(1)*dt1dr2+shpj(2)*dt2dr2+shpj(3)*dt3dr2)
|
||||
do m = 1,3
|
||||
drhodrm1(m) = 0.d0
|
||||
drhodrm2(m) = 0.d0
|
||||
drhodrm1(m) = dGamma2(i)*
|
||||
$ (t1i*drho1drm1(m)
|
||||
$ + t2i*drho2drm1(m)
|
||||
$ + t3i*drho3drm1(m))
|
||||
drhodrm2(m) = dGamma2(j)*
|
||||
$ (t1j*drho1drm2(m)
|
||||
$ + t2j*drho2drm2(m)
|
||||
$ + t3j*drho3drm2(m))
|
||||
enddo
|
||||
|
||||
c Compute derivatives wrt sij, but only if necessary
|
||||
if (dscrfcn(jn).ne.0.d0) then
|
||||
drho0ds1 = rhoa0j
|
||||
drho0ds2 = rhoa0i
|
||||
a1 = 2.d0/rij
|
||||
drho1ds1 = a1*rhoa1j*arg1i1
|
||||
drho1ds2 = a1*rhoa1i*arg1j1
|
||||
a2 = 2.d0/rij2
|
||||
drho2ds1 = a2*rhoa2j*arg1i2
|
||||
$ - 2.d0/3.d0*Arho2b(i)*rhoa2j
|
||||
drho2ds2 = a2*rhoa2i*arg1j2
|
||||
$ - 2.d0/3.d0*Arho2b(j)*rhoa2i
|
||||
a3 = 2.d0/rij3
|
||||
a3a = 6.d0/(5.d0*rij)
|
||||
drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3
|
||||
drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3
|
||||
|
||||
if (ialloy.eq.1) then
|
||||
|
||||
a1i = 0.d0
|
||||
a1j = 0.d0
|
||||
a2i = 0.d0
|
||||
a2j = 0.d0
|
||||
a3i = 0.d0
|
||||
a3j = 0.d0
|
||||
if ( tsq_ave(1,i) .ne. 0.d0 ) then
|
||||
a1i = rhoa0j/tsq_ave(1,i)
|
||||
endif
|
||||
if ( tsq_ave(1,j) .ne. 0.d0 ) then
|
||||
a1j = rhoa0i/tsq_ave(1,j)
|
||||
endif
|
||||
if ( tsq_ave(2,i) .ne. 0.d0 ) then
|
||||
a2i = rhoa0j/tsq_ave(2,i)
|
||||
endif
|
||||
if ( tsq_ave(2,j) .ne. 0.d0 ) then
|
||||
a2j = rhoa0i/tsq_ave(2,j)
|
||||
endif
|
||||
if ( tsq_ave(3,i) .ne. 0.d0 ) then
|
||||
a3i = rhoa0j/tsq_ave(3,i)
|
||||
endif
|
||||
if ( tsq_ave(3,j) .ne. 0.d0 ) then
|
||||
a3j = rhoa0i/tsq_ave(3,j)
|
||||
endif
|
||||
|
||||
dt1ds1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2)
|
||||
dt1ds2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2)
|
||||
dt2ds1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2)
|
||||
dt2ds2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2)
|
||||
dt3ds1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2)
|
||||
dt3ds2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2)
|
||||
|
||||
else if (ialloy.eq.2) then
|
||||
|
||||
dt1ds1 = 0.d0
|
||||
dt1ds2 = 0.d0
|
||||
dt2ds1 = 0.d0
|
||||
dt2ds2 = 0.d0
|
||||
dt3ds1 = 0.d0
|
||||
dt3ds2 = 0.d0
|
||||
|
||||
else
|
||||
|
||||
ai = 0.d0
|
||||
if( rho0(i) .ne. 0.d0 ) then
|
||||
ai = rhoa0j/rho0(i)
|
||||
end if
|
||||
aj = 0.d0
|
||||
if( rho0(j) .ne. 0.d0 ) then
|
||||
aj = rhoa0i/rho0(j)
|
||||
end if
|
||||
|
||||
dt1ds1 = ai*(t1_meam(eltj)-t1i)
|
||||
dt1ds2 = aj*(t1_meam(elti)-t1j)
|
||||
dt2ds1 = ai*(t2_meam(eltj)-t2i)
|
||||
dt2ds2 = aj*(t2_meam(elti)-t2j)
|
||||
dt3ds1 = ai*(t3_meam(eltj)-t3i)
|
||||
dt3ds2 = aj*(t3_meam(elti)-t3j)
|
||||
|
||||
endif
|
||||
|
||||
drhods1 = dGamma1(i)*drho0ds1
|
||||
$ + dGamma2(i)*
|
||||
$ (dt1ds1*rho1(i)+t1i*drho1ds1
|
||||
$ + dt2ds1*rho2(i)+t2i*drho2ds1
|
||||
$ + dt3ds1*rho3(i)+t3i*drho3ds1)
|
||||
$ - dGamma3(i)*
|
||||
$ (shpi(1)*dt1ds1+shpi(2)*dt2ds1+shpi(3)*dt3ds1)
|
||||
drhods2 = dGamma1(j)*drho0ds2
|
||||
$ + dGamma2(j)*
|
||||
$ (dt1ds2*rho1(j)+t1j*drho1ds2
|
||||
$ + dt2ds2*rho2(j)+t2j*drho2ds2
|
||||
$ + dt3ds2*rho3(j)+t3j*drho3ds2)
|
||||
$ - dGamma3(j)*
|
||||
$ (shpj(1)*dt1ds2+shpj(2)*dt2ds2+shpj(3)*dt3ds2)
|
||||
endif
|
||||
|
||||
c Compute derivatives of energy wrt rij, sij and rij(3)
|
||||
dUdrij = phip*sij
|
||||
$ + fp(i)*drhodr1 + fp(j)*drhodr2
|
||||
dUdsij = 0.d0
|
||||
if (dscrfcn(jn).ne.0.d0) then
|
||||
dUdsij = phi
|
||||
$ + fp(i)*drhods1 + fp(j)*drhods2
|
||||
endif
|
||||
do m = 1,3
|
||||
dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m)
|
||||
enddo
|
||||
|
||||
c Add the part of the force due to dUdrij and dUdsij
|
||||
|
||||
force = dUdrij*recip + dUdsij*dscrfcn(jn)
|
||||
do m = 1,3
|
||||
forcem = delij(m)*force + dUdrijm(m)
|
||||
f(m,i) = f(m,i) + forcem
|
||||
f(m,j) = f(m,j) - forcem
|
||||
enddo
|
||||
|
||||
c Tabulate per-atom virial as symmetrized stress tensor
|
||||
|
||||
if (vflag_atom.ne.0) then
|
||||
fi(1) = delij(1)*force + dUdrijm(1)
|
||||
fi(2) = delij(2)*force + dUdrijm(2)
|
||||
fi(3) = delij(3)*force + dUdrijm(3)
|
||||
v(1) = -0.5 * (delij(1) * fi(1))
|
||||
v(2) = -0.5 * (delij(2) * fi(2))
|
||||
v(3) = -0.5 * (delij(3) * fi(3))
|
||||
v(4) = -0.25 * (delij(1)*fi(2) + delij(2)*fi(1))
|
||||
v(5) = -0.25 * (delij(1)*fi(3) + delij(3)*fi(1))
|
||||
v(6) = -0.25 * (delij(2)*fi(3) + delij(3)*fi(2))
|
||||
|
||||
vatom(1,i) = vatom(1,i) + v(1)
|
||||
vatom(2,i) = vatom(2,i) + v(2)
|
||||
vatom(3,i) = vatom(3,i) + v(3)
|
||||
vatom(4,i) = vatom(4,i) + v(4)
|
||||
vatom(5,i) = vatom(5,i) + v(5)
|
||||
vatom(6,i) = vatom(6,i) + v(6)
|
||||
vatom(1,j) = vatom(1,j) + v(1)
|
||||
vatom(2,j) = vatom(2,j) + v(2)
|
||||
vatom(3,j) = vatom(3,j) + v(3)
|
||||
vatom(4,j) = vatom(4,j) + v(4)
|
||||
vatom(5,j) = vatom(5,j) + v(5)
|
||||
vatom(6,j) = vatom(6,j) + v(6)
|
||||
endif
|
||||
|
||||
c Now compute forces on other atoms k due to change in sij
|
||||
|
||||
if (sij.eq.0.d0.or.sij.eq.1.d0) goto 100
|
||||
do kn = 1,numneigh_full
|
||||
k = firstneigh_full(kn)
|
||||
eltk = fmap(type(k))
|
||||
if (k.ne.j.and.eltk.gt.0) then
|
||||
call dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2,
|
||||
$ ntype,type,fmap,x,scrfcn,fcpair)
|
||||
if (dsij1.ne.0.d0.or.dsij2.ne.0.d0) then
|
||||
force1 = dUdsij*dsij1
|
||||
force2 = dUdsij*dsij2
|
||||
do m = 1,3
|
||||
delik(m) = x(m,k) - x(m,i)
|
||||
deljk(m) = x(m,k) - x(m,j)
|
||||
enddo
|
||||
do m = 1,3
|
||||
f(m,i) = f(m,i) + force1*delik(m)
|
||||
f(m,j) = f(m,j) + force2*deljk(m)
|
||||
f(m,k) = f(m,k) - force1*delik(m)
|
||||
$ - force2*deljk(m)
|
||||
enddo
|
||||
|
||||
c Tabulate per-atom virial as symmetrized stress tensor
|
||||
|
||||
if (vflag_atom.ne.0) then
|
||||
fi(1) = force1*delik(1)
|
||||
fi(2) = force1*delik(2)
|
||||
fi(3) = force1*delik(3)
|
||||
fj(1) = force2*deljk(1)
|
||||
fj(2) = force2*deljk(2)
|
||||
fj(3) = force2*deljk(3)
|
||||
v(1) = -third * (delik(1)*fi(1) + deljk(1)*fj(1))
|
||||
v(2) = -third * (delik(2)*fi(2) + deljk(2)*fj(2))
|
||||
v(3) = -third * (delik(3)*fi(3) + deljk(3)*fj(3))
|
||||
v(4) = -sixth * (delik(1)*fi(2) + deljk(1)*fj(2) +
|
||||
$ delik(2)*fi(1) + deljk(2)*fj(1))
|
||||
v(5) = -sixth * (delik(1)*fi(3) + deljk(1)*fj(3) +
|
||||
$ delik(3)*fi(1) + deljk(3)*fj(1))
|
||||
v(6) = -sixth * (delik(2)*fi(3) + deljk(2)*fj(3) +
|
||||
$ delik(3)*fi(2) + deljk(3)*fj(2))
|
||||
|
||||
vatom(1,i) = vatom(1,i) + v(1)
|
||||
vatom(2,i) = vatom(2,i) + v(2)
|
||||
vatom(3,i) = vatom(3,i) + v(3)
|
||||
vatom(4,i) = vatom(4,i) + v(4)
|
||||
vatom(5,i) = vatom(5,i) + v(5)
|
||||
vatom(6,i) = vatom(6,i) + v(6)
|
||||
vatom(1,j) = vatom(1,j) + v(1)
|
||||
vatom(2,j) = vatom(2,j) + v(2)
|
||||
vatom(3,j) = vatom(3,j) + v(3)
|
||||
vatom(4,j) = vatom(4,j) + v(4)
|
||||
vatom(5,j) = vatom(5,j) + v(5)
|
||||
vatom(6,j) = vatom(6,j) + v(6)
|
||||
vatom(1,k) = vatom(1,k) + v(1)
|
||||
vatom(2,k) = vatom(2,k) + v(2)
|
||||
vatom(3,k) = vatom(3,k) + v(3)
|
||||
vatom(4,k) = vatom(4,k) + v(4)
|
||||
vatom(5,k) = vatom(5,k) + v(5)
|
||||
vatom(6,k) = vatom(6,k) + v(6)
|
||||
endif
|
||||
|
||||
endif
|
||||
endif
|
||||
c end of k loop
|
||||
enddo
|
||||
endif
|
||||
100 continue
|
||||
endif
|
||||
c end of j loop
|
||||
enddo
|
||||
|
||||
c else if elti=0, this is not a meam atom
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,111 +0,0 @@
|
|||
c
|
||||
c declaration in pair_meam.h:
|
||||
c
|
||||
c void meam_setup_global(int *, int *, double *, int *, double *, double *,
|
||||
c double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, double *,
|
||||
c double *, double *, int *);
|
||||
c
|
||||
c call in pair_meam.cpp:
|
||||
c
|
||||
c meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
|
||||
c alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
|
||||
c
|
||||
c
|
||||
|
||||
subroutine meam_setup_global(nelt, lat, z, ielement, atwt, alpha,
|
||||
$ b0, b1, b2, b3, alat, esub, asub,
|
||||
$ t0, t1, t2, t3, rozero, ibar)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer nelt, lat, ielement, ibar
|
||||
real*8 z, atwt, alpha, b0, b1, b2, b3
|
||||
real*8 alat, esub, asub, t0, t1, t2, t3
|
||||
real*8 rozero
|
||||
|
||||
dimension lat(nelt), ielement(nelt), ibar(nelt)
|
||||
dimension z(nelt), atwt(nelt), alpha(nelt)
|
||||
dimension b0(nelt), b1(nelt), b2(nelt), b3(nelt)
|
||||
dimension alat(nelt), esub(nelt), asub(nelt)
|
||||
dimension t0(nelt), t1(nelt), t2(nelt), t3(nelt), rozero(nelt)
|
||||
|
||||
integer i
|
||||
real*8 tmplat(maxelt)
|
||||
|
||||
neltypes = nelt
|
||||
|
||||
do i = 1,nelt
|
||||
|
||||
if (lat(i).eq.0) then
|
||||
lattce_meam(i,i) = 'fcc'
|
||||
else if (lat(i).eq.1) then
|
||||
lattce_meam(i,i) = 'bcc'
|
||||
else if (lat(i).eq.2) then
|
||||
lattce_meam(i,i) = 'hcp'
|
||||
else if (lat(i).eq.3) then
|
||||
lattce_meam(i,i) = 'dim'
|
||||
else if (lat(i).eq.4) then
|
||||
lattce_meam(i,i) = 'dia'
|
||||
else
|
||||
c unknown
|
||||
endif
|
||||
|
||||
Z_meam(i) = z(i)
|
||||
ielt_meam(i) = ielement(i)
|
||||
alpha_meam(i,i) = alpha(i)
|
||||
beta0_meam(i) = b0(i)
|
||||
beta1_meam(i) = b1(i)
|
||||
beta2_meam(i) = b2(i)
|
||||
beta3_meam(i) = b3(i)
|
||||
tmplat(i) = alat(i)
|
||||
Ec_meam(i,i) = esub(i)
|
||||
A_meam(i) = asub(i)
|
||||
t0_meam(i) = t0(i)
|
||||
t1_meam(i) = t1(i)
|
||||
t2_meam(i) = t2(i)
|
||||
t3_meam(i) = t3(i)
|
||||
rho0_meam(i) = rozero(i)
|
||||
ibar_meam(i) = ibar(i)
|
||||
|
||||
if (lattce_meam(i,i).eq.'fcc') then
|
||||
re_meam(i,i) = tmplat(i)/sqrt(2.d0)
|
||||
elseif (lattce_meam(i,i).eq.'bcc') then
|
||||
re_meam(i,i) = tmplat(i)*sqrt(3.d0)/2.d0
|
||||
elseif (lattce_meam(i,i).eq.'hcp') then
|
||||
re_meam(i,i) = tmplat(i)
|
||||
elseif (lattce_meam(i,i).eq.'dim') then
|
||||
re_meam(i,i) = tmplat(i)
|
||||
elseif (lattce_meam(i,i).eq.'dia') then
|
||||
re_meam(i,i) = tmplat(i)*sqrt(3.d0)/4.d0
|
||||
else
|
||||
c error
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
c Set some defaults
|
||||
rc_meam = 4.0
|
||||
delr_meam = 0.1
|
||||
attrac_meam(:,:) = 0.0
|
||||
repuls_meam(:,:) = 0.0
|
||||
Cmax_meam(:,:,:) = 2.8
|
||||
Cmin_meam(:,:,:) = 2.0
|
||||
ebound_meam(:,:) = (2.8d0**2)/(4.d0*(2.8d0-1.d0))
|
||||
delta_meam(:,:) = 0.0
|
||||
nn2_meam(:,:) = 0
|
||||
zbl_meam(:,:) = 1
|
||||
gsmooth_factor = 99.0
|
||||
augt1 = 1
|
||||
ialloy = 0
|
||||
mix_ref_t = 0
|
||||
emb_lin_neg = 0
|
||||
bkgd_dyn = 0
|
||||
erose_form = 0
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
|
@ -1,204 +0,0 @@
|
|||
c
|
||||
c do a sanity check on index parameters
|
||||
subroutine meam_checkindex(num,lim,nidx,idx,ierr)
|
||||
implicit none
|
||||
integer i,num,lim,nidx,idx(3),ierr
|
||||
|
||||
ierr = 0
|
||||
if (nidx.lt.num) then
|
||||
ierr = 2
|
||||
return
|
||||
endif
|
||||
|
||||
do i=1,num
|
||||
if ((idx(i).lt.1).or.(idx(i).gt.lim)) then
|
||||
ierr = 3
|
||||
return
|
||||
endif
|
||||
enddo
|
||||
end
|
||||
|
||||
c
|
||||
c Declaration in pair_meam.h:
|
||||
c
|
||||
c void meam_setup_param(int *, double *, int *, int *, int *);
|
||||
c
|
||||
c Call in pair_meam.cpp
|
||||
c
|
||||
c meam_setup_param(&which,&value,&nindex,index,&errorflag);
|
||||
c
|
||||
c
|
||||
c
|
||||
c The "which" argument corresponds to the index of the "keyword" array
|
||||
c in pair_meam.cpp:
|
||||
c
|
||||
c 0 = Ec_meam
|
||||
c 1 = alpha_meam
|
||||
c 2 = rho0_meam
|
||||
c 3 = delta_meam
|
||||
c 4 = lattce_meam
|
||||
c 5 = attrac_meam
|
||||
c 6 = repuls_meam
|
||||
c 7 = nn2_meam
|
||||
c 8 = Cmin_meam
|
||||
c 9 = Cmax_meam
|
||||
c 10 = rc_meam
|
||||
c 11 = delr_meam
|
||||
c 12 = augt1
|
||||
c 13 = gsmooth_factor
|
||||
c 14 = re_meam
|
||||
c 15 = ialloy
|
||||
c 16 = mixture_ref_t
|
||||
c 17 = erose_form
|
||||
c 18 = zbl_meam
|
||||
c 19 = emb_lin_neg
|
||||
c 20 = bkgd_dyn
|
||||
|
||||
subroutine meam_setup_param(which, value, nindex,
|
||||
$ index, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer which, nindex, index(3), errorflag
|
||||
real*8 value
|
||||
integer i1, i2
|
||||
|
||||
errorflag = 0
|
||||
|
||||
c 0 = Ec_meam
|
||||
if (which.eq.0) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
Ec_meam(index(1),index(2)) = value
|
||||
|
||||
c 1 = alpha_meam
|
||||
else if (which.eq.1) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
alpha_meam(index(1),index(2)) = value
|
||||
|
||||
c 2 = rho0_meam
|
||||
else if (which.eq.2) then
|
||||
call meam_checkindex(1,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
rho0_meam(index(1)) = value
|
||||
|
||||
c 3 = delta_meam
|
||||
else if (which.eq.3) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
delta_meam(index(1),index(2)) = value
|
||||
|
||||
c 4 = lattce_meam
|
||||
else if (which.eq.4) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
if (value.eq.0) then
|
||||
lattce_meam(index(1),index(2)) = "fcc"
|
||||
else if (value.eq.1) then
|
||||
lattce_meam(index(1),index(2)) = "bcc"
|
||||
else if (value.eq.2) then
|
||||
lattce_meam(index(1),index(2)) = "hcp"
|
||||
else if (value.eq.3) then
|
||||
lattce_meam(index(1),index(2)) = "dim"
|
||||
else if (value.eq.4) then
|
||||
lattce_meam(index(1),index(2)) = "dia"
|
||||
else if (value.eq.5) then
|
||||
lattce_meam(index(1),index(2)) = 'b1'
|
||||
else if (value.eq.6) then
|
||||
lattce_meam(index(1),index(2)) = 'c11'
|
||||
else if (value.eq.7) then
|
||||
lattce_meam(index(1),index(2)) = 'l12'
|
||||
else if (value.eq.8) then
|
||||
lattce_meam(index(1),index(2)) = 'b2'
|
||||
endif
|
||||
|
||||
c 5 = attrac_meam
|
||||
else if (which.eq.5) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
attrac_meam(index(1),index(2)) = value
|
||||
|
||||
c 6 = repuls_meam
|
||||
else if (which.eq.6) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
repuls_meam(index(1),index(2)) = value
|
||||
|
||||
c 7 = nn2_meam
|
||||
else if (which.eq.7) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
i1 = min(index(1),index(2))
|
||||
i2 = max(index(1),index(2))
|
||||
nn2_meam(i1,i2) = value
|
||||
|
||||
c 8 = Cmin_meam
|
||||
else if (which.eq.8) then
|
||||
call meam_checkindex(3,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
Cmin_meam(index(1),index(2),index(3)) = value
|
||||
|
||||
c 9 = Cmax_meam
|
||||
else if (which.eq.9) then
|
||||
call meam_checkindex(3,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
Cmax_meam(index(1),index(2),index(3)) = value
|
||||
|
||||
c 10 = rc_meam
|
||||
else if (which.eq.10) then
|
||||
rc_meam = value
|
||||
|
||||
c 11 = delr_meam
|
||||
else if (which.eq.11) then
|
||||
delr_meam = value
|
||||
|
||||
c 12 = augt1
|
||||
else if (which.eq.12) then
|
||||
augt1 = value
|
||||
|
||||
c 13 = gsmooth
|
||||
else if (which.eq.13) then
|
||||
gsmooth_factor = value
|
||||
|
||||
c 14 = re_meam
|
||||
else if (which.eq.14) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
re_meam(index(1),index(2)) = value
|
||||
|
||||
c 15 = ialloy
|
||||
else if (which.eq.15) then
|
||||
ialloy = value
|
||||
|
||||
c 16 = mixture_ref_t
|
||||
else if (which.eq.16) then
|
||||
mix_ref_t = value
|
||||
|
||||
c 17 = erose_form
|
||||
else if (which.eq.17) then
|
||||
erose_form = value
|
||||
|
||||
c 18 = zbl_meam
|
||||
else if (which.eq.18) then
|
||||
call meam_checkindex(2,maxelt,nindex,index,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
i1 = min(index(1),index(2))
|
||||
i2 = max(index(1),index(2))
|
||||
zbl_meam(i1,i2) = value
|
||||
|
||||
c 19 = emb_lin_neg
|
||||
else if (which.eq.19) then
|
||||
emb_lin_neg = value
|
||||
|
||||
c 20 = bkgd_dyn
|
||||
else if (which.eq.20) then
|
||||
bkgd_dyn = value
|
||||
|
||||
else
|
||||
errorflag = 1
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
|
|
@ -1 +0,0 @@
|
|||
../Install.py
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * Fortran Library for Reactive Force Field *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * CONTRIBUTING AUTHORS: Hansohl Cho(MIT), Aidan Thompson(SNL) *
|
||||
# * and Greg Wagner(SNL) *
|
||||
# * CONTACT: hansohl@mit.edu, athompson@sandia.gov, gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.gfortran
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = reax_connect.F reax_inout.F reax_lammps.F reax_poten.F reax_reac.F reax_charges.F
|
||||
|
||||
HEADERFILES = reax_defs.h *.blk
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libreax.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = g95
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F $(HEADERFILES)
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o $(LIB)
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * Fortran Library for Reactive Force Field *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * CONTRIBUTING AUTHORS: Hansohl Cho(MIT), Aidan Thompson(SNL) *
|
||||
# * and Greg Wagner(SNL) *
|
||||
# * CONTACT: hansohl@mit.edu, athompson@sandia.gov, gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.gfortran
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = reax_connect.F reax_inout.F reax_lammps.F reax_poten.F reax_reac.F reax_charges.F
|
||||
|
||||
HEADERFILES = reax_defs.h *.blk
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libreax.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = gfortran
|
||||
F90FLAGS = -O3 -fPIC -fno-second-underscore
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F $(HEADERFILES)
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o $(LIB)
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * Fortran Library for Reactive Force Field *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * CONTRIBUTING AUTHORS: Hansohl Cho(MIT), Aidan Thompson(SNL) *
|
||||
# * and Greg Wagner(SNL) *
|
||||
# * CONTACT: hansohl@mit.edu, athompson@sandia.gov, gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.ifort
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = reax_connect.F reax_inout.F reax_lammps.F reax_poten.F reax_reac.F reax_charges.F
|
||||
|
||||
HEADERFILES = reax_defs.h *.blk
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libreax.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = ifort
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F $(HEADERFILES)
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o $(LIB)
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
# Settings that the LAMMPS build will import when this package library is used
|
||||
|
||||
reax_SYSINC =
|
||||
reax_SYSLIB =
|
||||
reax_SYSPATH =
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
# Settings that the LAMMPS build will import when this package library is used
|
||||
|
||||
reax_SYSINC =
|
||||
reax_SYSLIB = -lgfortran
|
||||
reax_SYSPATH =
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
# Settings that the LAMMPS build will import when this package library is used
|
||||
|
||||
reax_SYSINC =
|
||||
reax_SYSLIB = -lifcore
|
||||
reax_SYSPATH =
|
||||
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * Fortran Library for Reactive Force Field *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * CONTRIBUTING AUTHORS: Hansohl Cho(MIT), Aidan Thompson(SNL) *
|
||||
# * and Greg Wagner(SNL) *
|
||||
# * CONTACT: hansohl@mit.edu, athompson@sandia.gov, gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.empty
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = reax_connect.F reax_inout.F reax_lammps.F reax_poten.F reax_reac.F reax_charges.F
|
||||
|
||||
HEADERFILES = reax_defs.h *.blk
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libreax.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = mpifort
|
||||
F90FLAGS = -O3 -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F $(HEADERFILES)
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o $(LIB)
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * Fortran Library for Reactive Force Field *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * CONTRIBUTING AUTHORS: Hansohl Cho(MIT), Aidan Thompson(SNL) *
|
||||
# * and Greg Wagner(SNL) *
|
||||
# * CONTACT: hansohl@mit.edu, athompson@sandia.gov, gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.pgf90
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = reax_connect.F reax_inout.F reax_lammps.F reax_poten.F reax_reac.F reax_charges.F
|
||||
|
||||
HEADERFILES = reax_defs.h *.blk
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libreax.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = pgf90
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F $(HEADERFILES)
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o $(LIB)
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * Fortran Library for Reactive Force Field *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * CONTRIBUTING AUTHORS: Hansohl Cho(MIT), Aidan Thompson(SNL) *
|
||||
# * and Greg Wagner(SNL) *
|
||||
# * CONTACT: hansohl@mit.edu, athompson@sandia.gov, gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.ifort
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = reax_connect.F reax_inout.F reax_lammps.F reax_poten.F reax_reac.F reax_charges.F
|
||||
|
||||
HEADERFILES = reax_defs.h *.blk
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libreax.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = mpif90
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F $(HEADERFILES)
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o $(LIB)
|
||||
|
|
@ -1 +0,0 @@
|
|||
Makefile.gfortran
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * Fortran Library for Reactive Force Field *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * CONTRIBUTING AUTHORS: Hansohl Cho(MIT), Aidan Thompson(SNL) *
|
||||
# * and Greg Wagner(SNL) *
|
||||
# * CONTACT: hansohl@mit.edu, athompson@sandia.gov, gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# which file will be copied to Makefile.lammps
|
||||
|
||||
EXTRAMAKE = Makefile.lammps.ifort
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = reax_connect.F reax_inout.F reax_lammps.F reax_poten.F reax_reac.F reax_charges.F
|
||||
|
||||
HEADERFILES = reax_defs.h *.blk
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libreax.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = mpif90
|
||||
F90FLAGS = -O -fPIC
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
@cp $(EXTRAMAKE) Makefile.lammps
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F $(HEADERFILES)
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o $(LIB)
|
||||
|
|
@ -1,78 +0,0 @@
|
|||
ReaxFF library
|
||||
|
||||
Aidan Thompson, Sandia National Labs
|
||||
athomps at sandia.gov
|
||||
Jan 2008
|
||||
|
||||
This library is an implementation of the ReaxFF potential,
|
||||
specifically designed to work with LAMMPS. It is derived from Adri van
|
||||
Duin's original serial code, with intervening incarnations in CMDF and
|
||||
GRASP.
|
||||
|
||||
-------------------------------------------------
|
||||
|
||||
This directory has source files to build a library that LAMMPS
|
||||
links against when using the REAX package.
|
||||
|
||||
This library must be built with a F90 compiler, before LAMMPS is
|
||||
built, so LAMMPS can link against it.
|
||||
|
||||
You can type "make lib-reax" from the src directory to see help on how
|
||||
to build this library via make commands, or you can do the same thing
|
||||
by typing "python Install.py" from within this directory, or you can
|
||||
do it manually by following the instructions below.
|
||||
|
||||
Build the library using one of the provided Makefile.* files or create
|
||||
your own, specific to your compiler and system. For example:
|
||||
|
||||
make -f Makefile.gfortran
|
||||
|
||||
When you are done building this library, two files should
|
||||
exist in this directory:
|
||||
|
||||
libreax.a the library LAMMPS will link against
|
||||
Makefile.lammps settings the LAMMPS Makefile will import
|
||||
|
||||
Makefile.lammps is created by the make command, by copying one of the
|
||||
Makefile.lammps.* files. See the EXTRAMAKE setting at the top of the
|
||||
Makefile.* files.
|
||||
|
||||
IMPORTANT: You must examine the final Makefile.lammps to insure it is
|
||||
correct for your system, else the LAMMPS build will likely fail.
|
||||
|
||||
Makefile.lammps has settings for 3 variables:
|
||||
|
||||
user-reax_SYSINC = leave blank for this package
|
||||
user-reax_SYSLIB = auxiliary F90 libs needed to link a F90 lib with
|
||||
a C++ program (LAMMPS) via a C++ compiler
|
||||
user-reax_SYSPATH = path(s) to where those libraries are
|
||||
|
||||
Because you have a F90 compiler on your system, you should have these
|
||||
libraries. But you will have to figure out which ones are needed and
|
||||
where they are. Examples of common configurations are in the
|
||||
Makefile.lammps.* files.
|
||||
|
||||
-------------------------------------------------
|
||||
|
||||
Additional build notes:
|
||||
|
||||
The include file reax_defs.h is used by both the ReaxFF library source
|
||||
files and the LAMMPS pair_reax.cpp source file (in package src/REAX).
|
||||
It contains dimensions of statically-allocated arrays created by the
|
||||
ReaxFF library. The size of these arrays must be set small enough to
|
||||
avoid exceeding the available machine memory, and large enough to fit
|
||||
the actual data generated by ReaxFF. If you change the values in
|
||||
reax_defs.h, you must first rebuild the library and then rebuild
|
||||
LAMMPS.
|
||||
|
||||
This library is called by functions in pair_reax.cpp. The C++ to
|
||||
FORTRAN function calls in pair_reax.cpp assume that FORTRAN object
|
||||
names are converted to C object names by appending an underscore
|
||||
character. This is generally the case, but on machines that do not
|
||||
conform to this convention, you will need to modify either the C++
|
||||
code or your compiler settings. The name conversion is handled by the
|
||||
preprocessor macro called FORTRAN in the file pair_reax_fortran.h,
|
||||
which is included by pair_reax.cpp. Different definitions of this
|
||||
macro can be obtained by adding a machine-specific macro definition to
|
||||
the CCFLAGS variable in your your LAMMPS Makefile e.g. -D_IBM. See
|
||||
pair_reax_fortran.h for more info.
|
||||
|
|
@ -1,116 +0,0 @@
|
|||
#include "reax_defs.h"
|
||||
implicit real*8 (a-h,o-z),integer(i-n)
|
||||
parameter (nneighmax=NNEIGHMAXDEF)
|
||||
parameter (nat=NATDEF) !Max number of atoms
|
||||
parameter (nattot=NATTOTDEF) !Max number of global atoms
|
||||
parameter (nsort=NSORTDEF) !Max number of atom types
|
||||
parameter (mbond=MBONDDEF) !Max number of bonds connected to one atom
|
||||
parameter (na1mx3=3*nat) !3*max number of atoms
|
||||
parameter (navib=NAVIBDEF) !for 2nd derivatives
|
||||
parameter (nbotym=NBOTYMDEF) !Max number of bond types
|
||||
parameter (nvatym=NVATYMDEF) !Max number of valency angle types
|
||||
parameter (ntotym=NTOTYMDEF) !Max number of torsion angle types
|
||||
parameter (nhbtym=NHBTYMDEF) !Max number of hydrogen bond types
|
||||
parameter (nodmtym=NODMTYMDEF) !Max number of off-diagonal Morse types
|
||||
parameter (nboallmax=NBOALLMAXDEF) !Max number of all bonds
|
||||
parameter (nbomax=NBOMAXDEF) !Max number of bonds
|
||||
parameter (nhbmax=NHBMAXDEF) !Max number of hydrogen bonds
|
||||
parameter (nvamax=NVAMAXDEF) !Max number of valency angles
|
||||
parameter (nopmax=NOPMAXDEF) !Max number of out of plane angles
|
||||
parameter (ntomax=NTOMAXDEF) !Max number of torsion angles
|
||||
parameter (npamax=NPAMAXDEF) !Max number of general parameters in force field
|
||||
parameter (nmolmax=NMOLMAXDEF) !Max number of molecules in system
|
||||
parameter (nmolset=NMOLSETDEF) !Max number of molecules in training set
|
||||
parameter (mrestra=MRESTRADEF) !Max number of restraints
|
||||
parameter (mtreg=MTREGDEF) !Max number of temperature regimes
|
||||
parameter (mtzone=MTZONEDEF) !Max number of temperature zones
|
||||
parameter (mvreg=MVREGDEF) !Max number of volume regimes
|
||||
parameter (mvzone=MVZONEDEF) !Max number of volume zones
|
||||
parameter (mereg=MEREGDEF) !Max number of electric field regimes
|
||||
parameter (mezone=MEZONEDEF) !Max number of electric field zones
|
||||
character*1 qr,qrset,qresi2
|
||||
character*2 qaset,qadd
|
||||
character*3 qresi1
|
||||
character*5 qlabel,qffty,qbgfaxes,qbgfsgn,qresi3
|
||||
character*20 qkeyw
|
||||
character*25 qfile
|
||||
character*40 qffield,qformat,qstrana2
|
||||
character*60 qremark,qremset,qmolset
|
||||
character*200 qstrana1
|
||||
common
|
||||
$/cbka/ dhbdc(nhbmax,3,3),cp(nat,3),
|
||||
$ cadd(nat,3),d2(3*navib,3*navib),
|
||||
$ veladd(3,nat),
|
||||
$ aold(3,nat),dic(3,nat),pvdw1(nsort,nsort),
|
||||
$ pvdw2(nsort,nsort),angimp(nat,6),
|
||||
$ yt(na1mx3),pt(na1mx3),gi(na1mx3),enmolset(nmolset),
|
||||
$ ai(na1mx3),bi(na1mx3),yi(na1mx3),pn(na1mx3),tbo(nat),
|
||||
$ chgbgf(nattot),
|
||||
$ abo2(nat),bor4(nat),bos(nbomax),
|
||||
$ eldef(nat),vradic(nat),
|
||||
$ vmo2(nat),
|
||||
$ ro(nbomax),dbondr(nbomax),
|
||||
$ dbosidr(nbomax),thgo(nopmax),
|
||||
$ elmol(nmolmax),
|
||||
$ elaf(nsort),vpq(nsort),
|
||||
$ rvdw(nsort),alf(nsort),eps(nsort),chat(nsort),
|
||||
$ rcore(nsort,nsort),ecore(nsort,nsort),acore(nsort,nsort),
|
||||
$ vlp2(nsort),
|
||||
$ valp2(nsort),vincr(nsort),
|
||||
$ vval3(nsort),
|
||||
$ vuncor(nbotym),
|
||||
$ vop(nsort),
|
||||
$ sigqeq(nsort),
|
||||
$ rrcha(mrestra),
|
||||
$ rmstra3(mrestra),
|
||||
$ rmcha(mrestra),
|
||||
$ rtcha(mrestra),rvcha(mrestra),
|
||||
$ v2bo(ntotym),v3bo(ntotym),
|
||||
$ eel,fctor,elr,
|
||||
$ presx2,presy2,presz2,
|
||||
$ tset2,
|
||||
$ enmol,formol,vvol,tpnrad,
|
||||
$ delvib,
|
||||
$ taut2,tincr,xmasmd,
|
||||
$ gdicmax,parc1,parc2,sumelec,
|
||||
$ xinh,fsnh,vqnh,snh,ham,errnh,sumhe,
|
||||
$ swa,swb2,swc0,swc1,swc2,swc3,swc4,swc5,swc6,
|
||||
$ swc7,plr,endpo2,ccpar,
|
||||
$ c4,estrmin,endpo,accincr,
|
||||
$ endpoold,xadd,yadd,zadd,addist,taddmol,
|
||||
$ Hug_E0, Hug_P0, Hug_V0, xImpVcm, shock_vel,
|
||||
$ shock_z_sep
|
||||
common
|
||||
$/cbka/
|
||||
$ ioop(nopmax,9),ifreqset(nmolset),
|
||||
$ ijk(nat,4),icgeopt(nmolset),
|
||||
$ irap(50),irdo(50,2),
|
||||
$ ityadd(nat),
|
||||
$ nmoloo(nat),iradic(nat),idef(nsort),nasort(nsort),
|
||||
$ ibgr1(nattot),ibgr2(nattot),idupc(6),
|
||||
$ imolsta(nat),
|
||||
$ ncent2(nbomax),irads,nrdd,nrddf,nbiolab,nuge,
|
||||
$ nbon2,npar,nodmty,ngnh,irac,nincrop,
|
||||
$ nboty,mdstep,
|
||||
$ nreac,
|
||||
$ nbonop,icelo2,
|
||||
$ iaddfreq,iveladd,invt,
|
||||
$ noop,ndtau,
|
||||
$ nelc3,nfc,nsav2,nmmax,ibh2,
|
||||
$ nmmaxold,nfcold,icellold,imodfile,
|
||||
$ icelo2old,inmov1,inmov2,nchaold,naa,nadattempt,
|
||||
$ nequi,iadj,
|
||||
$ ntest,nmm,
|
||||
$ nmolo5o,nradcount,nmollset,iflga,
|
||||
$ iperiod,ibgfversion,iremark,iconne,
|
||||
$ kx,ky,kz,iexco,iruid,ibity,nvlist,
|
||||
$ ityrad,iredo,iexx,iexy,iexz,ncellopt,
|
||||
$ ndata2,nprob,nit,i5758,ingeo,nmoloold,itemp,
|
||||
$ icgeo,ishock_type,isymm,
|
||||
$ qadd(nat),qlabel(nattot),qffty(nattot),qresi1(nattot),
|
||||
$ qresi2(nattot),qresi3(nattot),
|
||||
$ qremark(20),qformat(20),qr,qffield,
|
||||
$ qstrana1,qstrana2,qmolset(nmolset)
|
||||
***********************************************************************
|
||||
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
common
|
||||
$/cbkabo/ abo(nat)
|
||||
|
||||
|
||||
|
|
@ -1,3 +0,0 @@
|
|||
common
|
||||
$/cbkatomcoord/ id(nat,3),xmasat(nat),vel(3,nat),accel(3,nat)
|
||||
|
||||
|
|
@ -1,3 +0,0 @@
|
|||
common
|
||||
$/cbkbo/ bo(nbomax)
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
common
|
||||
$/cbkboncor/ dbosindc(nbomax,3,2*mbond+2),dbosidc(nbomax,3,2),
|
||||
$ bo131(nsort),bo132(nsort),bo133(nsort),
|
||||
$ ovc(nbotym),v13cor(nbotym)
|
||||
|
||||
|
|
@ -1,3 +0,0 @@
|
|||
common
|
||||
$/cbkbopi/ bopi(nbomax)
|
||||
|
||||
|
|
@ -1,3 +0,0 @@
|
|||
common
|
||||
$/cbkbopi2/ bopi2(nbomax)
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
common
|
||||
$/cbkbosi/ bosi(nbomax)
|
||||
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
common
|
||||
$/cbkc/ c(nat,3),cglobal(nattot,3),itag(nat),
|
||||
$chgglobal(nattot)
|
||||
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
common
|
||||
$/cbkch/ ch(nat)
|
||||
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
common
|
||||
$/cbkcha/ ech,syscha,chisys
|
||||
$ vfieldx,vfieldy,vfieldz,nmcharge,ioldchg
|
||||
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
common
|
||||
$/cbkcharmol/ vmcha(nmolmax),
|
||||
$ iat1mc(nmolmax),iat2mc(nmolmax)
|
||||
|
||||
|
|
@ -1,3 +0,0 @@
|
|||
common
|
||||
$/cbkchb/ chi(nsort),eta(nsort),gam(nsort)
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
common
|
||||
$/cbkconst/ dgrrdn,one,half,three,zero,caljou,rgasc,xjouca
|
||||
$ convmd
|
||||
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbkcovbon/ de2(nbotym),de3(nbotym),psi(nbotym),
|
||||
$ psp(nbotym),
|
||||
$ ltripstaball
|
||||
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
integer Lvirial,Latomvirial
|
||||
|
||||
common
|
||||
$/cbkd/ d(3,nat),estrain(nat)
|
||||
|
||||
common
|
||||
$/cbkvirial/ atomvirial(6,nat),virial(6),Lvirial,Latomvirial
|
||||
|
|
@ -1,3 +0,0 @@
|
|||
common
|
||||
$/cbkdbodc/ dbodc(nbomax,3,2)
|
||||
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
|
||||
|
||||
|
||||
common
|
||||
$/cbkdbopi2ndc/ dbopi2ndc(nbomax,3,2*mbond+2)
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/dbopidc/ dbopi2dc(nbomax,3,2),dbopidc(nbomax,3,2)
|
||||
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
|
||||
|
||||
|
||||
common
|
||||
$/dbopindc/ dbopindc(nbomax,3,2*mbond+2)
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbkdcell/ dcell(3,nat,27)
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbkdhdc/ dhdc(nvamax,3,3)
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
common
|
||||
$/cbkdistan/ axis(3),aaxh,baxh,caxh,iortho
|
||||
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbkdrdc/ drdc(nbomax,3,2)
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
|
||||
common
|
||||
$/cbkefield/ efix,efiy,efiz,c1
|
||||
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
|
||||
common
|
||||
$/cbkenergies/ eb,eoop,epen,estrc,deda(3),pressu,
|
||||
$ efi,elp,emol,ea,eres,et,eradbo,
|
||||
$ ev,eco,ecoa,ehb,sw,ew,ep,ekin
|
||||
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
character*5 qetype
|
||||
common
|
||||
$/cbkeregime/ qetype(mereg,mezone),nnereg(mereg),nerc,
|
||||
$ ereg(mereg,mezone),nitec(mereg)
|
||||
|
||||
|
|
@ -1,9 +0,0 @@
|
|||
character*2 qas
|
||||
common
|
||||
$/cbkff/ gamcco(nsort,nsort),vpar(npamax),vovun(nsort),
|
||||
$ stlp(nsort),aval(nsort),vlp1(nsort),
|
||||
$ vover(nbotym),valp1(nsort),
|
||||
$ vka(nvatym),qas(nsort),amas(nsort),e1(nbotym),
|
||||
$ valf(nsort),de1(nbotym),swb,nvs(nvatym,3),nso,nvaty
|
||||
|
||||
|
||||
|
|
@ -1,8 +0,0 @@
|
|||
|
||||
|
||||
|
||||
|
||||
common
|
||||
$/cbkfftorang/ v4(ntotym),vconj(ntotym),
|
||||
$ v1(ntotym),v2(ntotym),v3(ntotym)
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbkh/ h(nvamax)
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
common
|
||||
|
||||
$/cbkhbond/ hhb(nhbmax)
|
||||
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
|
||||
|
||||
|
||||
common
|
||||
$/cbkia/ ia(nat,mbond+3),iag(nat,mbond+3)
|
||||
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
|
||||
|
||||
|
||||
common
|
||||
$/cbkidbo/ idbo(nbomax,2*mbond+2),dbondc(nbomax,3,2*mbond+2),
|
||||
$ idbo1(nbomax)
|
||||
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbkimove/ imove(nattot)
|
||||
|
||||
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
|
||||
character*40 qruid
|
||||
common
|
||||
$/cbkinit/ tsetor,nzero,none,ntwo,nthree,qruid,systime,
|
||||
$ ustime,two,pi,avognr,axiss(3),pset,rdndgr
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbklonpar/ vlp(nat),dvlpdsbo(nat)
|
||||
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
|
||||
|
||||
|
||||
common
|
||||
$/cbkmolec/ nmolat2(nmolmax,nat),elmol2(nmolmax)
|
||||
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
|
||||
|
||||
|
||||
common
|
||||
$/cbknmolat/ nmolat(nmolmax,nat)
|
||||
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbknonbon/ gamwco(nsort,nsort),sw1,p3co(nsort,nsort),
|
||||
$ p2co(nsort,nsort),p1co(nsort,nsort)
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbknubon2/ nubon1(nat,mbond), nubon2(nat,mbond)
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
common
|
||||
$/cbknvlbo/ nvlbo(nneighmax*nat)
|
||||
|
||||
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
common
|
||||
$/cbknvlown/ nvlown(nneighmax*nat)
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
common
|
||||
$/cbkpairs/ nvl1(nneighmax*nat),nvl2(nneighmax*nat),nvpair,nvlself
|
||||
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
|
||||
common
|
||||
$/cbkpres/ presx,presy,presz
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
character*2 qa
|
||||
common
|
||||
$/cbkqa/ qa(nattot)
|
||||
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbkrbo/ rbo(nbomax),ibsym(nbomax),ib(nbomax,3)
|
||||
|
||||
|
|
@ -1,12 +0,0 @@
|
|||
common
|
||||
$/cbkrestr/ vkrv(mrestra),vrstra(mrestra),vkr2v(mrestra),
|
||||
$ dismacen(mrestra),rmstra1(mrestra),
|
||||
$ rrstra(mrestra),vkrst2(mrestra),
|
||||
$ rmstra2(mrestra),rmstrax(mrestra),rmstray(mrestra),
|
||||
$ rmstraz(mrestra),cmo(nat,3),vmo1(nat),trstra(mrestra),
|
||||
$ vkrt(mrestra),vkr2t(mrestra),vkrstr(mrestra),
|
||||
$ irstrav(mrestra,3),irstra(mrestra,2),itend(mrestra),
|
||||
$ irstram(mrestra,5),itstart(mrestra),irstrat(mrestra,4),
|
||||
$ imorph
|
||||
|
||||
|
||||
|
|
@ -1,13 +0,0 @@
|
|||
character*60 qmol
|
||||
common
|
||||
$/cbksrtbon1/ dbodr(nbomax),dbopidr(nbomax),
|
||||
$ dbopi2dr(nbomax),
|
||||
$ rob1(nsort,nsort),rob2(nsort,nsort),
|
||||
$ rob3(nsort,nsort),
|
||||
$ rat(nsort),rapt(nsort),vnq(nsort),bom(nbotym),
|
||||
$ pdp(nbotym),ptp(nbotym),pdo(nbotym),
|
||||
$ popi(nbotym),bop1(nbotym),bop2(nbotym),cutoff,
|
||||
$ nbs(nbotym,2),
|
||||
$ nsbma2,nsbmax,nboty2,nbonall,qfile(nmolset),qmol
|
||||
|
||||
|
||||
|
|
@ -1,10 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbksrthb/ vhb1(nhbtym),vhb2(nhbtym),rhb(nhbtym),
|
||||
$ dehb(nhbtym),ihb(nhbmax,8),nhb,
|
||||
$ nphb(nsort),nhbs(nhbtym,3),nhbty,hbcut,
|
||||
$ lhbnew
|
||||
|
||||
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbktorang/ dargtdc(ntomax,3,4),thg(ntomax)
|
||||
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
|
||||
|
||||
|
||||
common
|
||||
$/cbktorsion/ nts(ntotym,4),ntoty,ntor,it(ntomax,11)
|
||||
|
||||
|
||||
|
|
@ -1,8 +0,0 @@
|
|||
common
|
||||
$/cbktregime/ dttreg(mtreg,mtzone),tdamptreg(mtreg,mtzone),
|
||||
$ ia1treg(mtreg,mtzone),ia2treg(mtreg,mtzone),
|
||||
$ tsettreg(mtreg,mtzone),nntreg(mtreg),ittc(mtreg),
|
||||
$ nittc(mtreg),ifieldz,ifieldx,ifieldy,
|
||||
$ nrestra,nrestram,nrestrat,nrestrav,ntrc
|
||||
|
||||
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
|
||||
|
||||
common
|
||||
$/cbkvalence/ nval,iv(nvamax,6)
|
||||
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
character*5 qvtype
|
||||
common
|
||||
$/cbkvregime/ ivsca(mvreg,mvzone),dvvreg(mvreg,mvzone),
|
||||
$ nnvreg(mvreg),invrc,nitvc(mvreg),
|
||||
$ qvtype(mvreg,mvzone)
|
||||
|
||||
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
|
||||
common
|
||||
$/cellcoord/ tm11,tm21,tm31,tm22,tm32,tm33,angle(3),angles(3)
|
||||
|
||||
|
|
@ -1,12 +0,0 @@
|
|||
common
|
||||
$/control/ vrange,cutof2,cutof3,vlbora,tstep,range,taut,volcha,
|
||||
$ axis1,axis2,axis3,taup,
|
||||
$ icpres,nmethod,noutpt,inpt,napp,ianaly,ncha2,
|
||||
$ nrand,ntscale,itstep,ndebug,icentr,itrout,
|
||||
$ nchaudixmolo,itrans,nsav,
|
||||
$ nrep1,ncontrol,nhop2,nsav3,ngeofor,ifreq,
|
||||
$ nprevrun,maxstp,nvel,nsurp,ncons,
|
||||
$ ncha,icell,imolde,nchaud
|
||||
*********************************************************************************
|
||||
|
||||
|
||||
|
|
@ -1,23 +0,0 @@
|
|||
parameter (maxdat=5000)
|
||||
parameter (maxmdat=2500)
|
||||
parameter (maxkop=2500)
|
||||
character*80 qff
|
||||
character*60 qmdat
|
||||
character*100 qdatid
|
||||
character*2 qas2
|
||||
common
|
||||
$/opt/ fpar(7,nvatym,40),datopt(maxdat),caldat(maxdat),
|
||||
$ compdat(maxdat),weightdat(maxdat),
|
||||
$ devi(maxdat),vkop(maxkop),devkop(maxkop),sdy(3),
|
||||
$ valpar,valnew,change,vchange,
|
||||
$ molin(maxmdat,nsort),
|
||||
$ iboo(nbotym,2),idmo(nodmtym,2),ivao(nvatym,5),
|
||||
$ itoo(ntotym,7),ihbo(nhbtym,5),iheada(maxmdat),
|
||||
$ ndatm(maxmdat),iheada2(maxmdat),ichn(3),
|
||||
$ ikop1(maxkop),ikop2(maxkop),idat(maxdat),mu1(maxkop),
|
||||
$ mu2(maxkop),
|
||||
$ ndata,imam,iopt,iheatf,nkop,iagain,
|
||||
$ qdatid(maxdat),qmdat(maxmdat),qff(250),qas2(nsort)
|
||||
***********************************************************************
|
||||
|
||||
|
||||
|
|
@ -1,85 +0,0 @@
|
|||
**********************************************************************
|
||||
* *
|
||||
* REAXFF Reactive force field program *
|
||||
* *
|
||||
* Developed and written by Adri van Duin, duin@wag.caltech.edu *
|
||||
* *
|
||||
* Copyright (c) 2001-2010 California Institute of Technology *
|
||||
* *
|
||||
* This is an open-source program. Feel free to modify its *
|
||||
* contents. Please keep me informed of any useful modification *
|
||||
* or addition that you made. Please do not distribute this *
|
||||
* program to others; if people are interested in obtaining *
|
||||
* a copy of this program let them contact me first. *
|
||||
* *
|
||||
**********************************************************************
|
||||
************************************************************************
|
||||
|
||||
subroutine taper(r,r2)
|
||||
|
||||
************************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbkconst.blk"
|
||||
#include "cbkenergies.blk"
|
||||
#include "cbkinit.blk"
|
||||
#include "cbknonbon.blk"
|
||||
************************************************************************
|
||||
* *
|
||||
* Taper function for Coulomb interaction *
|
||||
* *
|
||||
************************************************************************
|
||||
r3=r2*r
|
||||
SW=SWC7*R3*R3*R+SWC6*R3*R3+SWC5*R3*R2+SWC4*R2*R2+SWC3*R3+SWC2*R2+
|
||||
$SWC1*R+SWC0
|
||||
SW1=7.0D0*SWC7*R3*R3+6.0D0*SWC6*R3*R2+5.0D0*SWC5*R2*R2+
|
||||
$4.0D0*SWC4*R3+THREE*SWC3*R2+TWO*SWC2*R+SWC1
|
||||
return
|
||||
end
|
||||
************************************************************************
|
||||
************************************************************************
|
||||
|
||||
subroutine tap7th
|
||||
|
||||
************************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbkc.blk"
|
||||
#include "cbkconst.blk"
|
||||
#include "cbkd.blk"
|
||||
#include "cbkff.blk"
|
||||
#include "cbkh.blk"
|
||||
#include "control.blk"
|
||||
|
||||
************************************************************************
|
||||
* *
|
||||
* 7th order taper function setup *
|
||||
* *
|
||||
************************************************************************
|
||||
c$$$ if (ndebug.eq.1) then
|
||||
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
||||
c$$$ write (65,*) 'In tap7th'
|
||||
c$$$ call timer(65)
|
||||
c$$$ close (65)
|
||||
c$$$ end if
|
||||
D1=SWB-SWA
|
||||
D7=D1**7.0D0
|
||||
SWA2=SWA*SWA
|
||||
SWA3=SWA2*SWA
|
||||
SWB2=SWB*SWB
|
||||
SWB3=SWB2*SWB
|
||||
|
||||
************************************************************************
|
||||
* 7th order taper function *
|
||||
************************************************************************
|
||||
|
||||
SWC7= 20.0D0/D7
|
||||
SWC6= -70.0D0*(SWA+SWB)/D7
|
||||
SWC5= 84.0D0*(SWA2+3.0D0*SWA*SWB+SWB2)/D7
|
||||
SWC4= -35.0D0*(SWA3+9.0D0*SWA2*SWB+9.0D0*SWA*SWB2+SWB3)/D7
|
||||
SWC3= 140.0D0*(SWA3*SWB+3.0D0*SWA2*SWB2+SWA*SWB3)/D7
|
||||
SWC2=-210.0D0*(SWA3*SWB2+SWA2*SWB3)/D7
|
||||
SWC1= 140.0D0*SWA3*SWB3/D7
|
||||
SWC0=(-35.0D0*SWA3*SWB2*SWB2+21.0D0*SWA2*SWB3*SWB2-
|
||||
$7.0D0*SWA*SWB3*SWB3+SWB3*SWB3*SWB)/D7
|
||||
|
||||
return
|
||||
END
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,70 +0,0 @@
|
|||
#define PORTABLECOMMENTFLAG
|
||||
#ifndef PORTABLECOMMENTFLAG
|
||||
// This is just a way to have portable comments
|
||||
// for both C++ and FORTRAN preprocessing.
|
||||
/* ///:EOH~ */
|
||||
/* */
|
||||
/* This file contains array dimension parameters for all the main */
|
||||
/* ReaxFF data structures, some of which need to be directly accessed */
|
||||
/* by Grasp C++ functions. If they are set too small, the calculation */
|
||||
/* will run out of allocated memory. If they are set too big, the machine */
|
||||
/* will not be able to allocate enough memory. */
|
||||
/* */
|
||||
|
||||
/* NNEIGHMAXDEF = Max number of neighbors / NATDEF */
|
||||
/* NATDEF = Max number of atoms */
|
||||
/* NATTOTDEF = Max number of global atoms */
|
||||
/* NSORTDEF = Max number of atom types */
|
||||
/* MBONDDEF = Max number of bonds connected to one atom */
|
||||
/* NAVIBDEF = for 2nd derivatives */
|
||||
/* NBOTYMDEF = Max number of bond types */
|
||||
/* NVATYMDEF = Max number of valency angle types */
|
||||
/* NTOTYMDEF = Max number of torsion angle types */
|
||||
/* NHBTYMDEF = Max number of hydrogen bond types */
|
||||
/* NODMTYMDEF = Max number of off-diagonal Morse types */
|
||||
/* NBOALLMAXDEF = Max number of all bonds */
|
||||
/* NBOMAXDEF = Max number of bonds */
|
||||
/* NHBMAXDEF = Max number of hydrogen bonds */
|
||||
/* NVAMAXDEF = Max number of valency angles */
|
||||
/* NOPMAXDEF = Max number of out of plane angles */
|
||||
/* NTOMAXDEF = Max number of torsion angles */
|
||||
/* NPAMAXDEF = Max number of general parameters in force field */
|
||||
/* NMOLMAXDEF = Max number of molecules in system */
|
||||
/* NMOLSETDEF = Max number of molecules in training set */
|
||||
/* MRESTRADEF = Max number of restraints */
|
||||
/* MTREGDEF = Max number of temperature regimes */
|
||||
/* MTZONEDEF = Max number of temperature zones */
|
||||
/* MVREGDEF = Max number of volume regimes */
|
||||
/* MVZONEDEF = Max number of volume zones */
|
||||
/* MEREGDEF = Max number of electric field regimes */
|
||||
/* MEZONEDEF = Max number of electric field zones */
|
||||
#endif
|
||||
|
||||
#define NNEIGHMAXDEF 120
|
||||
#define NATDEF 40000
|
||||
#define NATTOTDEF 39744
|
||||
#define NSORTDEF 20
|
||||
#define MBONDDEF 20
|
||||
#define NAVIBDEF 50
|
||||
#define NBOTYMDEF 200
|
||||
#define NVATYMDEF 200
|
||||
#define NTOTYMDEF 200
|
||||
#define NHBTYMDEF 200
|
||||
#define NODMTYMDEF 20
|
||||
#define NBOALLMAXDEF 180000
|
||||
#define NBOMAXDEF 90000
|
||||
#define NHBMAXDEF 400000
|
||||
#define NVAMAXDEF 300000
|
||||
#define NOPMAXDEF 00010
|
||||
#define NTOMAXDEF 65000
|
||||
#define NPAMAXDEF 50
|
||||
#define NMOLMAXDEF 2000
|
||||
#define NMOLSETDEF 1500
|
||||
#define MRESTRADEF 100
|
||||
#define MTREGDEF 100
|
||||
#define MTZONEDEF 5
|
||||
#define MVREGDEF 100
|
||||
#define MVZONEDEF 6
|
||||
#define MEREGDEF 100
|
||||
#define MEZONEDEF 3
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,392 +0,0 @@
|
|||
**********************************************************************
|
||||
* *
|
||||
* REAXFF Reactive force field program *
|
||||
* *
|
||||
* Developed and written by Adri van Duin, duin@wag.caltech.edu *
|
||||
* *
|
||||
* Copyright (c) 2001-2010 California Institute of Technology *
|
||||
* *
|
||||
* This is an open-source program. Feel free to modify its *
|
||||
* contents. Please keep me informed of any useful modification *
|
||||
* or addition that you made. Please do not distribute this *
|
||||
* program to others; if people are interested in obtaining *
|
||||
* a copy of this program let them contact me first. *
|
||||
* *
|
||||
**********************************************************************
|
||||
**********************************************************************
|
||||
|
||||
subroutine getswb(swb_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbkff.blk"
|
||||
real*8 swb_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of swb *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
swb_tmp = swb
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getswa(swa_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbkff.blk"
|
||||
real*8 swa_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of swa *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
swa_tmp = swa
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getvrange(vrange_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "control.blk"
|
||||
real*8 vrange_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of vrange *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
vrange_tmp = vrange
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnvlist(nvlist_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
integer nvlist_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nvlist *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nvlist_tmp = nvlist
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getvlbora(vlbora_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "control.blk"
|
||||
real*8 vlbora_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of vlbora *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
vlbora_tmp = vlbora
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnval(nval_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbkvalence.blk"
|
||||
integer nval_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nval *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nval_tmp = nval
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getntor(ntor_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbktorsion.blk"
|
||||
integer ntor_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of ntor *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
ntor_tmp = ntor
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnhb(nhb_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbksrthb.blk"
|
||||
integer nhb_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nhb *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nhb_tmp = nhb
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnbonall(nbonall_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbksrtbon1.blk"
|
||||
integer nbonall_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nbonall *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nbonall_tmp = nbonall
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnneighmax(nneighmax_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
integer nneighmax_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nneighmax *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nneighmax_tmp = nneighmax
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnat(nat_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
integer nat_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nat *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nat_tmp = nat
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnattot(nattot_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
integer nattot_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nattot *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nattot_tmp = nattot
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnsort(nsort_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
integer nsort_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nsort *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nsort_tmp = nsort
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getmbond(mbond_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
integer mbond_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of mbond *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
mbond_tmp = mbond
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnso(nso_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbkff.blk"
|
||||
integer nso_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nso *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nso_tmp = nso
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine setngeofor(ngeofor_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "control.blk"
|
||||
integer ngeofor_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Set value of ngeofor
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
ngeofor = ngeofor_tmp
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnsbmax(nsbmax_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbksrtbon1.blk"
|
||||
integer nsbmax_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nsbmax *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nsbmax_tmp = nsbmax
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getnsbma2(nsbma2_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "cbksrtbon1.blk"
|
||||
integer nsbma2_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of nsbma2 *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
nsbma2_tmp = nsbma2
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
**********************************************************************
|
||||
|
||||
subroutine getcutof3(cutof3_tmp)
|
||||
|
||||
**********************************************************************
|
||||
#include "cbka.blk"
|
||||
#include "control.blk"
|
||||
real*8 cutof3_tmp
|
||||
|
||||
**********************************************************************
|
||||
* *
|
||||
* Report the value of cutof3 *
|
||||
* *
|
||||
**********************************************************************
|
||||
|
||||
cutof3_tmp = cutof3
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
Some files were not shown because too many files have changed in this diff Show More
Loading…
Reference in New Issue