ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

more LATTE refactoring and debugging

Browse Source
This commit is contained in:
Steve Plimpton
2023-02-02 18:16:45 -07:00
parent 12079c32de
commit 751e33c70e
12 changed files with 367 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
examples/QUANTUM/LATTE/README.new
View File

@ -103,12 +103,54 @@ Step 4: run the water AIMD problem for a few steps
% lmp_mpi -mdi "-name LMP -role DRIVER -method TCP -port 8021" -log log.water.latte.aimd.tcp.1 -in in.water.latte.aimd &
% python latte_mdi.py -mdi "-name LATTE -role ENGINE -method TCP -port 8021 -hostname localhost"
% python latte_mdi.py -mdi "-name LATTE -role ENGINE -method TCP -port 8021 -hostname localhost" latte.in.water
# Run with MPI: 1 proc each
% mpirun -np 1 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" -log log.water.latte.aimd.mpi.1 -in in.water.latte.aimd : -np 1 python latte_mdi.py -mdi "-name LATTE -role ENGINE -method MPI"
% mpirun -np 1 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" -log log.water.latte.aimd.mpi.1 -in in.water.latte.aimd : -np 1 python latte_mdi.py -mdi "-name LATTE -role ENGINE -method MPI" latte.in.water
# Run in plugin mode: 1 proc
% lmp_mpi -mdi "-name LMP -role DRIVER -method LINK -plugin_path /home/sjplimp/lammps/git/examples/QUANTUM/LATTE" -log log.water.latte.aimd.plugin.1 -in in.water.latte.aimd.plugin
---------------------------------
---------------------------------
Step 5: run the UO2 AIMD problem for a few steps
% cd ~/lammps/examples/QUANTUM/LATTE
# Run with TCP: 1 proc each
% lmp_mpi -mdi "-name LMP -role DRIVER -method TCP -port 8021" -log log.uo2.latte.aimd.tcp.1 -in in.uo2.latte.aimd &
% python latte_mdi.py -mdi "-name LATTE -role ENGINE -method TCP -port 8021 -hostname localhost" latte.in.uo2
# Run with MPI: 1 proc each
% mpirun -np 1 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" -log log.uo2.latte.aimd.mpi.1 -in in.uo2.latte.aimd : -np 1 python latte_mdi.py -mdi "-name LATTE -role ENGINE -method MPI" latte.in.uo2
# Run in plugin mode: 1 proc
% lmp_mpi -mdi "-name LMP -role DRIVER -method LINK -plugin_path /home/sjplimp/lammps/git/examples/QUANTUM/LATTE" -log log.uo2.latte.aimd.plugin.1 -in in.uo2.latte.aimd.plugin
---------------------------------
---------------------------------
Step 6: run the UO2 series problem
% cd ~/lammps/examples/QUANTUM/LATTE
# Run with TCP: 1 proc each
% lmp_mpi -mdi "-name LMP -role DRIVER -method TCP -port 8021" -log log.uo2.latte.series.tcp.1 -in in.uo2.latte.series &
% python latte_mdi.py -mdi "-name LATTE -role ENGINE -method TCP -port 8021 -hostname localhost" latte.in.uo2
# Run with MPI: 1 proc each
% mpirun -np 1 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" -log log.uo2.latte.series.mpi.1 -in in.uo2.latte.series : -np 1 python latte_mdi.py -mdi "-name LATTE -role ENGINE -method MPI" latte.in.uo2
# Run in plugin mode: 1 proc
% lmp_mpi -mdi "-name LMP -role DRIVER -method LINK -plugin_path /home/sjplimp/lammps/git/examples/QUANTUM/LATTE" -log log.uo2.latte.series.plugin.1 -in in.uo2.latte.series.plugin

26
examples/QUANTUM/LATTE/in.uo2.latte.aimd Normal file
View File

@ -0,0 +1,26 @@
# AIMD test of two UO2 molecules with LATTE in MDI stand-alone mode
units metal
atom_style full
atom_modify sort 0 0.0
read_data 2uo2.lmp
velocity all create 300.0 87287 loop geom
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all mdi/qm virial yes elements U O
thermo_style custom step temp pe etotal press
thermo 1
#dump 1 all custom 1 dump.aimd.mpi &
# id type x y z vx vy vz fx fy fz
run 20

27
examples/QUANTUM/LATTE/in.uo2.latte.aimd.plugin Normal file
View File

@ -0,0 +1,27 @@
# AIMD test of two UO2 molecules with LATTE in MDI plugin mode
units metal
atom_style full
atom_modify sort 0 0.0
read_data 2uo2.lmp
velocity all create 300.0 87287 loop geom
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all mdi/qm virial yes elements U O
thermo_style custom step temp pe etotal press
thermo 1
#dump 1 all custom 1 dump.aimd.plugin &
# id type x y z vx vy vz fx fy fz
mdi plugin latte_mdi mdi "-role ENGINE -name LATTE -method LINK" &
extra latte.in.uo2 command "run 20"

37
examples/QUANTUM/LATTE/in.uo2.latte.series Normal file
View File

@ -0,0 +1,37 @@
# Series of single-point calcs of 2,3,4 UO2 molecules
# with LATTE in MDI stand-alone mode
variable files index 2uo2 3uo2 4uo2
mdi connect
label LOOP
units metal
atom_style full
atom_modify sort 0 0.0
read_data ${files}.lmp
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
fix 1 all mdi/qm virial yes elements U O connect no
thermo_style custom step temp pe etotal press
thermo 1
run 0
write_dump all custom dump.series.mpi.${files} &
id type x y z fx fy fz modify sort id
clear
next files
jump SELF LOOP
mdi exit

32
examples/QUANTUM/LATTE/in.uo2.latte.series.plugin Normal file
View File

@ -0,0 +1,32 @@
# Series of single-point calcs of 2,3,4 UO2 molecules
# with LATTE in MDI plugin mode
variable files index 2uo2 3uo2 4uo2
label LOOP
units metal
atom_style full
atom_modify sort 0 0.0
read_data ${files}.lmp
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
fix 1 all mdi/qm virial yes elements U O
thermo_style custom step temp pe etotal press
thermo 1
mdi plugin latte_mdi mdi "-role ENGINE -name LATTE -method LINK" &
extra latte.in.uo2 command "run 0"
write_dump all custom dump.series.plugin.${files} &
id type x y z fx fy fz modify sort id
clear
next files
jump SELF LOOP

6
examples/QUANTUM/LATTE/in.water.latte.aimd
View File

@ -29,7 +29,7 @@ timestep 0.00025
fix 1 all nve
fix 2 all mdi/qm virial yes elements 8 1
fix 2 all mdi/qm virial yes elements O H
#fix 2 all latte
fix_modify 2 energy yes
@ -37,5 +37,5 @@ thermo_style custom step temp pe etotal press
# dynamics
thermo 1
run 10
thermo 10
run 100

4
examples/QUANTUM/LATTE/in.water.latte.aimd.plugin
View File

@ -29,7 +29,7 @@ timestep 0.00025
fix 1 all nve
fix 2 all mdi/qm virial yes elements 8 1
fix 2 all mdi/qm virial yes elements O H
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
@ -39,4 +39,4 @@ thermo_style custom step temp pe etotal press
thermo 10
mdi plugin latte_mdi mdi "-role ENGINE -name LATTE -method LINK" &
command "run 100"
extra latte.in.water command "run 100"

186
examples/QUANTUM/LATTE/latte_mdi.py
View File

@ -1,15 +1,30 @@
# MDI wrapper on LATTE code
import sys,time
# native LATTE units are Angstroms and eV
import sys,os,time
from ctypes import *
import numpy as np
from mpi4py import MPI
import MDI_Library as mdi
# conversions of atomic number to element symbol
# --------------------------------------------
atomic_number_to_symbol = {1: 'H', 6: 'C', 7: 'N', 8: 'O', 17: 'Cl'}
ELEMENTS = [
'H' , 'He', 'Li', 'Be', 'B' , 'C' , 'N' , 'O' , 'F' , 'Ne',
'Na', 'Mg', 'Al', 'Si', 'P' , 'S' , 'Cl', 'Ar', 'K' , 'Ca',
'Sc', 'Ti', 'V' , 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn',
'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y' , 'Zr',
'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn',
'Sb', 'Te', 'I' , 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd',
'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb',
'Lu', 'Hf', 'Ta', 'W' , 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg',
'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th',
'Pa', 'U' , 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm',
'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds',
'Rg', 'Cn', 'Nh', 'Fl', 'Mc', 'Lv', 'Ts', 'Og',
]
# --------------------------------------------
# global data
@ -28,6 +43,7 @@ mode = AIMD
libname = "liblatte.so"
liblatte = None
infile = None
# QM inputs
@ -45,6 +61,10 @@ qm_natoms = 0
qm_elements = None
qm_coords = None
qm_potential = None
qm_velocity = None
qm_ntypes = 0
qm_types = None
qm_mass = None
mm_natoms = 0
mm_coords = None
@ -73,9 +93,16 @@ def error(txt):
# --------------------------------------------
def options(other_options):
if len(other_options) != 0:
error("No args currently used by LATTE wrapper")
global infile
if len(other_options) != 1:
error("No filename provided to LATTE wrapper")
infile = other_options[0]
# copy infile to latte.in as required by LATTE
os.system("cp %s latte.in" % infile)
# --------------------------------------------
# operate as an engine
# --------------------------------------------
@ -382,17 +409,21 @@ def receive_mm_charges(mdicomm):
# --------------------------------------------
def allocate(which):
global qm_types,qm_elements,qm_coords,qm_potential,qm_forces,qm_charges
global qm_elements,qm_coords,qm_potential,qm_forces,qm_charges
global qm_types,qm_mass,qm_velocity
global mm_coords,mm_charges,mm_forces
if which == "qm":
n = qm_natoms
qm_types = np.empty(n,dtype=np.int32)
qm_elements = np.empty(n,dtype=np.int32)
qm_coords = np.empty((n,3))
qm_potential = np.empty(n)
qm_forces = np.empty((n,3))
qm_charges = np.empty(n)
qm_types = np.empty(n,dtype=np.int32)
qm_velocity = np.empty((n,3))
qm_velocity.fill(0.0)
if which == "mm":
n = mm_natoms
@ -443,6 +474,7 @@ def evaluate():
error("QM atom properties not fully specified")
if mode == QMMM and not flag_qm_potential:
error("QM atom properties not fully specified")
set_types_masses()
# LATTE inputs
# box, qm_coords must be converted to Angstroms
@ -471,27 +503,14 @@ def evaluate():
#flags_latte[4] = vflag_atom and thermo_virial; # virial, 0 for no
#flags_latte[5] = neighflag; # 1 to pass neighbor list to LATTE, 0 for no
# hard-wired values for 8-water problem
qm_types = np.empty(qm_natoms,dtype=np.int32)
for i in range(0,qm_natoms,3): qm_types[i] = 1
for i in range(1,qm_natoms,3): qm_types[i] = 2
for i in range(2,qm_natoms,3): qm_types[i] = 2
qm_mass = [15.995, 1.008]
# setup ctypes args for liblatte.latte() function
flags_latte = 6*[0]
c_flags_latte = (c_int*6)(*flags_latte)
c_qm_natoms = c_int(qm_natoms)
qm_ntypes = 2
c_qm_ntypes = c_int(qm_ntypes)
c_qm_mass = (c_double*qm_ntypes)(*qm_mass)
boxlo = [0.0,0.0,0.0]
boxhi = [box_A[0],box_A[4],box_A[8]]
xy = box_A[3]
@ -509,10 +528,7 @@ def evaluate():
c_qm_pe = c_double(qm_pe)
qm_velocity = np.empty((qm_natoms,3))
qm_velocity.fill(0.0)
timestep = 0.00025
timestep = 0.0
c_timestep = c_double(timestep)
latte_stress = np.empty(6)
@ -520,11 +536,8 @@ def evaluate():
latte_error = 0
c_latte_error = c_bool(latte_error)
# QMMM with QM and MM atoms
# NOTE: need qm_velocity and timestep and mass and types ?
# all of these are addresses of scalars for Fortran ?
if mode == QMMM:
error("QMMM not yet supported with LATTE")
@ -535,53 +548,51 @@ def evaluate():
#print("Calling LATTE ...")
time1 = time.time()
#print("flags_latte",c_flags_latte[0:6])
#print("qm_natoms",c_qm_natoms.value)
#print("qm_coords",qm_coords_A)
#print("qm_types",qm_types)
#print("qm_ntypes",c_qm_ntypes.value)
#print("qm_mass",c_qm_mass[0:2])
#print("boxlo",c_boxlo[0:3])
#print("boxhi",c_boxhi[0:3])
#print("xy",c_xy.value)
#print("xz",c_xz.value)
#print("yz",c_yz.value)
#print("maxiter",c_maxiter.value)
#print("timestep",c_timestep.value)
#print("new_system",c_new_system.value)
#print("A",c_qm_natoms.value)
#print("B",qm_coords_A)
#print("C",qm_types)
#print("D",c_qm_ntypes.value)
#print("E",qm_mass)
#print("F",c_boxlo[0],c_boxlo[1],c_boxlo[2])
#print("G",c_boxhi[0],c_boxhi[1],c_boxhi[2])
#print("H",c_xy.value,c_xz.value,c_yz.value)
#print("I",c_maxiter.value)
#print("J",qm_velocity)
#print("K",c_timestep.value)
#print("L",c_new_system.value)
liblatte.\
latte(c_flags_latte,byref(c_qm_natoms),qm_coords_A,
qm_types,byref(c_qm_ntypes),c_qm_mass,
qm_types,byref(c_qm_ntypes),qm_mass,
c_boxlo,c_boxhi,byref(c_xy),byref(c_xz),byref(c_yz),
qm_forces,byref(c_maxiter),byref(c_qm_pe),
qm_velocity,byref(c_timestep),latte_stress,
byref(c_new_system),byref(c_latte_error))
# NOTE: check latte_error return?
latte_error = c_latte_error.value
time2 = time.time()
#print("LATTE time:",time2-time1,latte_error)
# LATTE computes extensive stress values
# MDI stress values are intensive
# so divide latte_stress by box volume in Angstroms
# more generally could use volume = (a x b) dot c for (a,b,c) edge vectors
qm_pe = c_qm_pe.value
latte_stress /= box_A[0]*box_A[4]*box_A[8]
#print("LATTE STRESS:",latte_stress)
time2 = time.time()
#print("DONE LATTE",latte_error,time2-time1)
# conversion of LATTE outputs
# conversion of LATTE outputs to MDI units
# qm_pe from eV to Hartrees
# qm_forces from eV/A to Hartrees/Bohr
# latte_stress from eV/A^3 to Hartrees/Bohr^3
# also need to divide stress by box volume to make it intensive
# NOTE: why box volume in Angs ? b/c that is what LATTE did ?
# did this above
ev_to_hartree = mdi.MDI_Conversion_factor("electron_volt","hartree")
angstrom_to_bohr = mdi.MDI_Conversion_factor("angstrom","bohr")
qm_pe *= ev_to_hartree
qm_forces *= ev_to_hartree / angstrom_to_bohr
latte_stress *= ev_to_hartree / (angstrom_to_bohr*angstrom_to_bohr*angstrom_to_bohr)
latte_stress *= ev_to_hartree / \
(angstrom_to_bohr*angstrom_to_bohr*angstrom_to_bohr)
qm_stress[0] = latte_stress[0]
qm_stress[4] = latte_stress[1]
@ -590,13 +601,6 @@ def evaluate():
qm_stress[2] = qm_stress[6] = latte_stress[4]
qm_stress[5] = qm_stress[7] = latte_stress[5]
#print("CONVERT:",ev_to_hartree / (angstrom_to_bohr*angstrom_to_bohr*angstrom_to_bohr))
#print("QM STRESS:",qm_stress)
#print("PE",qm_pe)
#print("FORCE",qm_forces)
#print("STRESS",qm_stress)
# clear flags for all MDI commands for next QM evaluation
flag_qm_natoms = flag_mm_natoms = 0
@ -625,12 +629,68 @@ def latte_load():
liblatte.latte.restype = None
liblatte.latte.argtypes = \
[POINTER(c_int), POINTER(c_int), nparray, npvector_int, POINTER(c_int), POINTER(c_double),
[POINTER(c_int), POINTER(c_int), nparray, npvector_int, POINTER(c_int), npvector_double,
POINTER(c_double), POINTER(c_double), POINTER(c_double),
POINTER(c_double), POINTER(c_double),
nparray, POINTER(c_int), POINTER(c_double), nparray,
POINTER(c_double), npvector_double, POINTER(c_int), POINTER(c_bool)]
# --------------------------------------------
# set ntypes, qm_types, qm_mass from qm_elements and electrons.dat file
# required inputs for LATTE latte() function
# --------------------------------------------
def set_types_masses():
global qm_ntypes,qm_types,qm_mass
# read latte.in file for PARAMPATH to electrons.dat file
lines = open(infile,'r').readlines()
for line in lines:
if "PARAMPATH=" not in line: continue
words = line.split()
parampath = words[1][1:-1]
if parampath[-1] == '/': efile = parampath + "electrons.dat"
else: efile = parampath + "/electrons.dat"
# extract symbol/mass pairings from electrons.dat file
# symbol2mass = dict for mass of each chemical symbol
lines = open(efile,'r').readlines()
symbol2mass = {}
for line in lines[2:]:
line = line.strip()
if not line: continue
words = line.split()
if words[0] in symbol2mass: continue
else: symbol2mass[words[0]] = float(words[7])
# symbol_list = list of unique chemical symbols
# qm_ntypes = # of unique symbols
symbol_list = []
for i in range(qm_natoms):
symbol = ELEMENTS[qm_elements[i]-1]
if symbol not in symbol2mass:
error("Atom element symbol not in LATTE electrons.dat file")
if symbol in symbol_list: continue
symbol_list.append(symbol)
qm_ntypes = len(symbol_list)
# qm_types = chemcial type of each atom (1 to qm_types)
# qm_mass = mass of each atom type
qm_mass = np.empty(qm_ntypes)
for i in range(qm_natoms):
symbol = ELEMENTS[qm_elements[i]-1]
itype = symbol_list.index(symbol)
qm_types[i] = itype + 1
qm_mass[itype] = symbol2mass[symbol]
# --------------------------------------------
# function called by MDI driver
# only when it invokes pyscf_mdi.py as a plugin

2
examples/QUANTUM/PySCF/pyscf_mdi.py
View File

@ -38,8 +38,6 @@ from pyscf.pbc.dft import RKS as RKS_pbc
atomic_number_to_radius = {1: 0.32, 6: 0.75, 7: 0.71, 8: 0.63, 17: 0.99}
# ELEMENTS is from pyscf/pyscf/data/elements.py
ELEMENTS = [
'H' , 'He', 'Li', 'Be', 'B' , 'C' , 'N' , 'O' , 'F' , 'Ne',
'Na', 'Mg', 'Al', 'Si', 'P' , 'S' , 'Cl', 'Ar', 'K' , 'Ca',

35
src/MDI/fix_mdi_qm.cpp
View File

@ -25,7 +25,7 @@ using namespace FixConst;
enum { NATIVE, REAL, METAL }; // LAMMPS units which MDI supports
#define MAXELEMENT 103 // used elsewhere in MDI package
#define MAXELEMENT 118
/* ---------------------------------------------------------------------- */
@ -87,14 +87,32 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
error->all(FLERR, "Illegal fix mdi/qm command");
iarg += 2;
} else if (strcmp(arg[iarg], "elements") == 0) {
const char *symbols[] = {
"H" , "He", "Li", "Be", "B" , "C" , "N" , "O" , "F" , "Ne",
"Na", "Mg", "Al", "Si", "P" , "S" , "Cl", "Ar", "K" , "Ca",
"Sc", "Ti", "V" , "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
"Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y" , "Zr",
"Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
"Sb", "Te", "I" , "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
"Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
"Lu", "Hf", "Ta", "W" , "Re", "Os", "Ir", "Pt", "Au", "Hg",
"Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
"Pa", "U" , "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
"Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds",
"Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og",
};
int ntypes = atom->ntypes;
if (iarg + ntypes + 1 > narg) error->all(FLERR, "Illegal fix mdi/qm command");
delete[] elements;
elements = new int[ntypes + 1];
for (int i = 1; i <= ntypes; i++) {
elements[i] = utils::inumeric(FLERR, arg[iarg + i], false, lmp);
if (elements[i] < 1 || elements[i] > MAXELEMENT)
error->all(FLERR, "Illegal fix mdi/qm command");
int anum;
for (anum = 0; anum < MAXELEMENT; anum++)
if (strcmp(arg[iarg + i],symbols[anum]) == 0) break;
if (anum == MAXELEMENT)
error->all(FLERR,"Invalid chemical element in fix mdi/qm command");
elements[i] = anum + 1;
}
iarg += ntypes + 1;
} else
@ -274,8 +292,8 @@ void FixMDIQM::init()
error->all(FLERR, "MDI: Engine has wrong atom count and does not support >NATOMS command");
}
int elements_exists;
int types_exists;
int elements_exists,types_exists;
ierr = MDI_Check_command_exists("@DEFAULT", ">ELEMENTS", mdicomm, &elements_exists);
if (ierr) error->all(FLERR, "MDI: >ELEMENTS command check");
MPI_Bcast(&elements_exists, 1, MPI_INT, 0, world);
@ -288,6 +306,7 @@ void FixMDIQM::init()
send_elements();
else if (types_exists)
send_types();
send_box();
}
@ -531,6 +550,7 @@ void FixMDIQM::send_elements()
/* ----------------------------------------------------------------------
send simulation box size and shape to MDI engine
only send CELL_DISPL if engine supports it
------------------------------------------------------------------------- */
void FixMDIQM::send_box()
@ -564,7 +584,8 @@ void FixMDIQM::send_box()
cell[7] = domain->yz;
cell[8] = domain->boxhi[2] - domain->boxlo[2];
// convert the cell units to bohr
// convert cell units to bohr
for (int icell = 0; icell < 9; icell++) cell[icell] *= lmp2mdi_length;
ierr = MDI_Send(cell, 9, MDI_DOUBLE, mdicomm);

31
src/MDI/fix_mdi_qmmm.cpp
View File

@ -31,7 +31,7 @@ using namespace FixConst;
enum { NATIVE, REAL, METAL }; // LAMMPS units which MDI supports
enum { DIRECT, POTENTIAL }; // mode of QMMM coupling
#define MAXELEMENT 103 // used elsewhere in MDI package
#define MAXELEMENT 118
// prototype for non-class compare function for sorting QM IDs
@ -43,9 +43,6 @@ FixMDIQMMM::FixMDIQMMM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
// check requirements for LAMMPS to work with MDI for a QMMM engine
if (domain->dimension == 2)
error->all(FLERR,"Fix mdi/qmmm requires 3d simulation");
if (!atom->tag_enable)
error->all(FLERR,"Fix mdi/qmmm requires atom IDs be defined");
@ -95,14 +92,32 @@ FixMDIQMMM::FixMDIQMMM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
error->all(FLERR, "Illegal fix mdi/qmmm command");
iarg += 2;
} else if (strcmp(arg[iarg], "elements") == 0) {
const char *symbols[] = {
"H" , "He", "Li", "Be", "B" , "C" , "N" , "O" , "F" , "Ne",
"Na", "Mg", "Al", "Si", "P" , "S" , "Cl", "Ar", "K" , "Ca",
"Sc", "Ti", "V" , "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
"Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y" , "Zr",
"Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
"Sb", "Te", "I" , "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
"Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
"Lu", "Hf", "Ta", "W" , "Re", "Os", "Ir", "Pt", "Au", "Hg",
"Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
"Pa", "U" , "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
"Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds",
"Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og",
};
int ntypes = atom->ntypes;
if (iarg + ntypes + 1 > narg) error->all(FLERR, "Illegal fix mdi/qmmm command");
delete[] elements;
elements = new int[ntypes + 1];
for (int i = 1; i <= ntypes; i++) {
elements[i] = utils::inumeric(FLERR, arg[iarg + i], false, lmp);
if (elements[i] < 1 || elements[i] > MAXELEMENT)
error->all(FLERR, "Illegal fix mdi/qmmm command");
for (int i = 1; i <= ntypes; i++) {
int anum;
for (anum = 0; anum < MAXELEMENT; anum++)
if (strcmp(arg[iarg + 1],symbols[anum]) == 0) break;
if (anum == MAXELEMENT)
error->all(FLERR,"Invalid chemical element in fix mdi/qmmm command");
elements[i] = anum + 1;
}
iarg += ntypes + 1;
} else

26
src/MDI/mdi_engine.cpp
View File

@ -54,7 +54,7 @@ enum { DEFAULT, MD, OPT }; // top-level MDI engine modes
enum { TYPE, CHARGE, MASS, COORD, VELOCITY, FORCE, ADDFORCE };
#define MAXELEMENT 103 // used elsewhere in MDI package
#define MAXELEMENT 118
/* ----------------------------------------------------------------------
trigger LAMMPS to start acting as an MDI engine
@ -81,14 +81,32 @@ MDIEngine::MDIEngine(LAMMPS *_lmp, int narg, char **arg) : Pointers(_lmp)
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg], "elements") == 0) {
const char *symbols[] = {
"H" , "He", "Li", "Be", "B" , "C" , "N" , "O" , "F" , "Ne",
"Na", "Mg", "Al", "Si", "P" , "S" , "Cl", "Ar", "K" , "Ca",
"Sc", "Ti", "V" , "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
"Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y" , "Zr",
"Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
"Sb", "Te", "I" , "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
"Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
"Lu", "Hf", "Ta", "W" , "Re", "Os", "Ir", "Pt", "Au", "Hg",
"Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
"Pa", "U" , "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
"Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds",
"Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og",
};
int ntypes = atom->ntypes;
delete[] elements;
elements = new int[ntypes + 1];
if (iarg + ntypes + 1 > narg) error->all(FLERR, "Illegal mdi engine command");
for (int i = 1; i <= ntypes; i++) {
elements[i] = utils::inumeric(FLERR, arg[iarg + i], false, lmp);
if (elements[i] < 0 || elements[i] > MAXELEMENT)
error->all(FLERR, "Illegal mdi engine command");
int anum;
for (anum = 0; anum < MAXELEMENT; anum++)
if (strcmp(arg[iarg + 1],symbols[anum]) == 0) break;
if (anum == MAXELEMENT)
error->all(FLERR,"Invalid chemical element in mdi engine command");
elements[i] = anum + 1;
}
iarg += ntypes + 1;
} else