add support for NWChem in examples/QUANTUM
This commit is contained in:
Steve Plimpton
119
examples/QUANTUM/NWChem/README
Normal file
119
examples/QUANTUM/NWChem/README
Normal file
@ -0,0 +1,119 @@
|
|||||||
|
# Test runs of QMMM with LAMMPS and NWChem
|
||||||
|
|
||||||
|
Step 1: build LAMMPS
|
||||||
|
Step 2: download/build the MDI code coupling package
|
||||||
|
Step 3: download/build NWChem PWDFT
|
||||||
|
Step 4: run 2-water QMMM problem for a few steps
|
||||||
|
|
||||||
|
---------------------------------
|
||||||
|
---------------------------------
|
||||||
|
|
||||||
|
Step 1: build LAMMPS
|
||||||
|
|
||||||
|
The molecule package is needed for the 2-water test
|
||||||
|
problem. Copy the final LAMMPS executable into the
|
||||||
|
examples/QUANTUM/NWChem directory.
|
||||||
|
|
||||||
|
Traditional make:
|
||||||
|
|
||||||
|
% cd ~/lammps/lib/mdi
|
||||||
|
% python Install.py -m mpi
|
||||||
|
% cd ~/lammps/src
|
||||||
|
% make yes-mdi yes-molecule
|
||||||
|
% make -j mpi
|
||||||
|
% cp lmp_mpi ~/lammps/examples/QUANTUM/NWChem
|
||||||
|
|
||||||
|
CMake:
|
||||||
|
|
||||||
|
% cd ~/lammps
|
||||||
|
% mkdir build; cd build
|
||||||
|
% cmake -D PKG_MDI=yes -D PKG_MOLECULE=yes ../cmake
|
||||||
|
% make -j
|
||||||
|
% cp lmp ~/lammps/examples/QUANTUM/NWChem/lmp_mpi
|
||||||
|
|
||||||
|
---------------------------------
|
||||||
|
---------------------------------
|
||||||
|
|
||||||
|
Step 2: install the MDI code coupling package
|
||||||
|
|
||||||
|
(a) grab the MDI Git repo
|
||||||
|
|
||||||
|
% mkdir mdi; cd mdi
|
||||||
|
% git clone git@github.com:MolSSI-MDI/MDI_Library.git git
|
||||||
|
|
||||||
|
(b) build MDI
|
||||||
|
|
||||||
|
% cd mdi/git
|
||||||
|
% mkdir build; cd build
|
||||||
|
% cmake ..
|
||||||
|
% make -j
|
||||||
|
|
||||||
|
(c) Add something similar to the following to your .bashrc or .cshrc
|
||||||
|
file so that Python can find MDI:
|
||||||
|
|
||||||
|
For bash:
|
||||||
|
|
||||||
|
% export PYTHONPATH="$PYTHONPATH:/home/sjplimp/mdi/git"
|
||||||
|
% hash -r
|
||||||
|
|
||||||
|
For (t)csh:
|
||||||
|
|
||||||
|
% setenv PYTHONPATH ${PYTHONPATH}:/home/sjplimp/mdi/git
|
||||||
|
% rehash
|
||||||
|
|
||||||
|
(d) Check that you can import the 3 Python modules which the script
|
||||||
|
that wraps PySCF will need:
|
||||||
|
|
||||||
|
% python
|
||||||
|
>>> import numpy as np
|
||||||
|
>>> from mpi4py import MPI
|
||||||
|
>>> import MDI_Library as mdi
|
||||||
|
|
||||||
|
---------------------------------
|
||||||
|
---------------------------------
|
||||||
|
|
||||||
|
Step 3: download/build NWChem PWDFT
|
||||||
|
|
||||||
|
(a) grab the PWDFT Git repo
|
||||||
|
|
||||||
|
% mkdir nwchem; cd nwchem
|
||||||
|
% git clone git@github.com:ebylaska/PWDFT.git PWDFT
|
||||||
|
|
||||||
|
(b) build PWDFT
|
||||||
|
|
||||||
|
% cd nwchem/PWDFT
|
||||||
|
% cd build_library; rm -r *
|
||||||
|
% cmake ../Nwpw -DMAKE_LIBRARY=true -DCMAKE_POSITION_INDEPENDENT_CODE=ON
|
||||||
|
% make -j # should produce libpwdft.so in build_library
|
||||||
|
|
||||||
|
(c) Add something similar to the following to your .bashrc or .cshrc
|
||||||
|
file so that the Python wrapper script can find libpwdft.so:
|
||||||
|
|
||||||
|
For bash:
|
||||||
|
|
||||||
|
% export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:/home/sjplimp/nwchem/PWDFT/build_library"
|
||||||
|
% hash -r
|
||||||
|
|
||||||
|
For (t)csh:
|
||||||
|
|
||||||
|
% setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:/home/sjplimp/nwchem/PWDFT/build_library
|
||||||
|
% rehash
|
||||||
|
|
||||||
|
---------------------------------
|
||||||
|
---------------------------------
|
||||||
|
|
||||||
|
Step 4: run the 2-water QMMM problem for a few steps
|
||||||
|
|
||||||
|
% cd ~/lammps/examples/QUANTUM/NWChem
|
||||||
|
|
||||||
|
lmp_mpi -mdi "-name LMP -role DRIVER -method TCP -port 8021" -log log.water.nwchem.qmmm.tcp.1 -in in.water.nwchem.qmmm &
|
||||||
|
|
||||||
|
python nwchem_mdi.py -mdi "-name NWChem -role ENGINE -method TCP -port 8021 -hostname localhost" template.water.nw water.dimer.nw log.water.pwdft.qmmm.tcp.1
|
||||||
|
|
||||||
|
# Run with MPI: 1 proc each
|
||||||
|
|
||||||
|
mpirun -np 1 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" -log log.water.nwchem.qmmm.mpi.1 -in in.water.nwchem.qmmm : -np 1 python nwchem_mdi.py -mdi "-name NWChem -role ENGINE -method MPI" template.water.nw water.dimer.nw log.water.pwdft.qmmm.mpi.1
|
||||||
|
|
||||||
|
# Run in plugin mode: 1 proc
|
||||||
|
|
||||||
|
lmp_mpi -mdi "-name LMP -role DRIVER -method LINK -plugin_path /home/sjplimp/work_qm" -log log.water.nwchem.qmmm.plugin.1 -in in.water.nwchem.qmmm.plugin
|
||||||
48
examples/QUANTUM/NWChem/data.water.nwchem
Normal file
48
examples/QUANTUM/NWChem/data.water.nwchem
Normal file
@ -0,0 +1,48 @@
|
|||||||
|
LAMMPS data file for water dimer in large box - for QMMM with different types
|
||||||
|
|
||||||
|
6 atoms
|
||||||
|
4 atom types
|
||||||
|
4 bonds
|
||||||
|
1 bond types
|
||||||
|
2 angles
|
||||||
|
1 angle types
|
||||||
|
|
||||||
|
-6.879301 6.879301 xlo xhi
|
||||||
|
-6.879301 6.879301 ylo yhi
|
||||||
|
-6.879301 6.879301 zlo zhi
|
||||||
|
|
||||||
|
Masses
|
||||||
|
|
||||||
|
1 15.99491
|
||||||
|
2 1.008
|
||||||
|
3 15.99491
|
||||||
|
4 1.008
|
||||||
|
|
||||||
|
Bond Coeffs
|
||||||
|
|
||||||
|
1 554.25 1.0
|
||||||
|
|
||||||
|
Angle Coeffs
|
||||||
|
|
||||||
|
1 47.744 109.4
|
||||||
|
|
||||||
|
Atoms
|
||||||
|
|
||||||
|
1 1 1 -0.8476 0.161560 -0.052912 0.033173
|
||||||
|
2 1 2 0.4238 0.803054 0.369132 -0.511660
|
||||||
|
3 1 2 0.4238 -0.325571 -0.669574 -0.488560
|
||||||
|
4 2 3 -0.8476 0.021259 0.506771 2.831278
|
||||||
|
5 2 4 0.4238 -0.721039 1.083100 2.758378
|
||||||
|
6 2 4 0.4238 0.158220 0.181883 1.945696
|
||||||
|
|
||||||
|
Bonds
|
||||||
|
|
||||||
|
1 1 1 2
|
||||||
|
2 1 1 3
|
||||||
|
3 1 4 5
|
||||||
|
4 1 4 6
|
||||||
|
|
||||||
|
Angles
|
||||||
|
|
||||||
|
1 1 2 1 3
|
||||||
|
2 1 5 4 6
|
||||||
70
examples/QUANTUM/NWChem/in.water.nwchem.qmmm
Normal file
70
examples/QUANTUM/NWChem/in.water.nwchem.qmmm
Normal file
@ -0,0 +1,70 @@
|
|||||||
|
# QMMM with NWChem
|
||||||
|
|
||||||
|
units real
|
||||||
|
atom_style full
|
||||||
|
|
||||||
|
bond_style harmonic
|
||||||
|
angle_style harmonic
|
||||||
|
|
||||||
|
read_data data.water.nwchem
|
||||||
|
|
||||||
|
# QM atoms are 1st water
|
||||||
|
# MM atoms are 2nd water
|
||||||
|
|
||||||
|
group qm molecule 1
|
||||||
|
group mm molecule 2
|
||||||
|
|
||||||
|
# remove bonds/angles between QM atoms
|
||||||
|
|
||||||
|
delete_bonds qm multi remove special
|
||||||
|
|
||||||
|
# pair style must define stand-alone short-range Coulombics
|
||||||
|
# must specify mixing explicitly b/c hybrid/overlay
|
||||||
|
# QM O,H = types 1,2
|
||||||
|
# MM O,H = types 3,4
|
||||||
|
# QM O,H atoms do not LJ interact with each other
|
||||||
|
# only MM O atoms LJ interact with other b/c MM H is zero
|
||||||
|
# MM/QM O do LJ interact with each other, same as pair of MM O atoms
|
||||||
|
# MM O and QM H do LJ interact with each other with non-zero H epsilon = 0.044
|
||||||
|
# geometric mixing for epsilon, arithmetic for sigma
|
||||||
|
# this is to provide stability for QM H atoms
|
||||||
|
|
||||||
|
# mixing only for MM-O/QM-O and MM-O/QM-H
|
||||||
|
|
||||||
|
pair_style hybrid/overlay lj/cut 6.0 coul/cut 6.0
|
||||||
|
pair_coeff 1 1 lj/cut 0.0 3.165558
|
||||||
|
pair_coeff 2 2 lj/cut 0.0 0.7
|
||||||
|
pair_coeff 3 3 lj/cut 0.155394 3.165558
|
||||||
|
pair_coeff 4 4 lj/cut 0.0 0.7
|
||||||
|
pair_coeff 1 3 lj/cut 0.155394 3.165558
|
||||||
|
pair_coeff 2 3 lj/cut 0.08268818537130924 1.932779
|
||||||
|
pair_coeff * * coul/cut
|
||||||
|
|
||||||
|
neighbor 1.0 bin
|
||||||
|
neigh_modify delay 0 every 1 check yes
|
||||||
|
|
||||||
|
# QMMM dynamics
|
||||||
|
|
||||||
|
timestep 0.1
|
||||||
|
|
||||||
|
fix 1 all nve
|
||||||
|
|
||||||
|
fix 2 qm mdi/qmmm potential elements 8 1 8 1
|
||||||
|
#fix 2 qm nwchem template.water.nw water.dimer2.nw &
|
||||||
|
# log.pwdft.water.dimer O H O H
|
||||||
|
fix_modify 2 energy yes
|
||||||
|
|
||||||
|
thermo_style custom step cpu temp ke evdwl ecoul epair emol elong &
|
||||||
|
f_2 pe etotal press
|
||||||
|
|
||||||
|
# convert dump file forces to Hartree/Bohr for comparison to NWChem
|
||||||
|
|
||||||
|
variable fx atom fx/1185.8
|
||||||
|
variable fy atom fy/1185.8
|
||||||
|
variable fz atom fz/1185.8
|
||||||
|
|
||||||
|
dump 1 all custom 1 dump.water.nwchem.qmmm id x y z q v_fx v_fy v_fz
|
||||||
|
dump_modify 1 sort id format float "%20.16g"
|
||||||
|
|
||||||
|
thermo 1
|
||||||
|
run 2
|
||||||
72
examples/QUANTUM/NWChem/in.water.nwchem.qmmm.plugin
Normal file
72
examples/QUANTUM/NWChem/in.water.nwchem.qmmm.plugin
Normal file
@ -0,0 +1,72 @@
|
|||||||
|
# QMMM with NWChem
|
||||||
|
|
||||||
|
units real
|
||||||
|
atom_style full
|
||||||
|
|
||||||
|
bond_style harmonic
|
||||||
|
angle_style harmonic
|
||||||
|
|
||||||
|
read_data data.water.nwchem
|
||||||
|
|
||||||
|
# QM atoms are 1st water
|
||||||
|
# MM atoms are 2nd water
|
||||||
|
|
||||||
|
group qm molecule 1
|
||||||
|
group mm molecule 2
|
||||||
|
|
||||||
|
# remove bonds/angles between QM atoms
|
||||||
|
|
||||||
|
delete_bonds qm multi remove special
|
||||||
|
|
||||||
|
# pair style must define stand-alone short-range Coulombics
|
||||||
|
# must specify mixing explicitly b/c hybrid/overlay
|
||||||
|
# QM O,H = types 1,2
|
||||||
|
# MM O,H = types 3,4
|
||||||
|
# QM O,H atoms do not LJ interact with each other
|
||||||
|
# only MM O atoms LJ interact with other b/c MM H is zero
|
||||||
|
# MM/QM O do LJ interact with each other, same as pair of MM O atoms
|
||||||
|
# MM O and QM H do LJ interact with each other with non-zero H epsilon = 0.044
|
||||||
|
# geometric mixing for epsilon, arithmetic for sigma
|
||||||
|
# this is to provide stability for QM H atoms
|
||||||
|
|
||||||
|
# mixing only for MM-O/QM-O and MM-O/QM-H
|
||||||
|
|
||||||
|
pair_style hybrid/overlay lj/cut 6.0 coul/cut 6.0
|
||||||
|
pair_coeff 1 1 lj/cut 0.0 3.165558
|
||||||
|
pair_coeff 2 2 lj/cut 0.0 0.7
|
||||||
|
pair_coeff 3 3 lj/cut 0.155394 3.165558
|
||||||
|
pair_coeff 4 4 lj/cut 0.0 0.7
|
||||||
|
pair_coeff 1 3 lj/cut 0.155394 3.165558
|
||||||
|
pair_coeff 2 3 lj/cut 0.08268818537130924 1.932779
|
||||||
|
pair_coeff * * coul/cut
|
||||||
|
|
||||||
|
neighbor 1.0 bin
|
||||||
|
neigh_modify delay 0 every 1 check yes
|
||||||
|
|
||||||
|
# QMMM dynamics
|
||||||
|
|
||||||
|
timestep 0.1
|
||||||
|
|
||||||
|
fix 1 all nve
|
||||||
|
|
||||||
|
fix 2 qm mdi/qmmm potential elements 8 1 8 1
|
||||||
|
fix_modify 2 energy yes
|
||||||
|
|
||||||
|
thermo_style custom step cpu temp ke evdwl ecoul epair emol elong &
|
||||||
|
f_2 pe etotal press
|
||||||
|
|
||||||
|
# convert dump file forces to Hartree/Bohr for comparison to NWChem
|
||||||
|
|
||||||
|
variable fx atom fx/1185.8
|
||||||
|
variable fy atom fy/1185.8
|
||||||
|
variable fz atom fz/1185.8
|
||||||
|
|
||||||
|
dump 1 all custom 1 dump.water.dimer.qmmm.plugin &
|
||||||
|
id x y z q v_fx v_fy v_fz
|
||||||
|
dump_modify 1 sort id format float "%20.16g"
|
||||||
|
|
||||||
|
thermo 1
|
||||||
|
|
||||||
|
mdi plugin nwchem_mdi mdi "-role ENGINE -name NWChem -method LINK" &
|
||||||
|
extra "template.water.nw water.dimer.nw log.water.pwdft.qmmm.plugin.1" &
|
||||||
|
command "run 2"
|
||||||
662
examples/QUANTUM/NWChem/nwchem_mdi.py
Normal file
662
examples/QUANTUM/NWChem/nwchem_mdi.py
Normal file
@ -0,0 +1,662 @@
|
|||||||
|
# MDI wrapper on NWChem PWDFT code
|
||||||
|
|
||||||
|
# NOTE: Qs or issues to still address
|
||||||
|
# test if works for both AIMD and QMMM
|
||||||
|
# can series of problem be run via lib interface input_filename() ?
|
||||||
|
# how does PBC vs non-PBC work, just box size in NWC input file
|
||||||
|
# or maybe other settings in that file?
|
||||||
|
# can NWChem return stress?
|
||||||
|
# can NWChem do DIRECT mode?
|
||||||
|
# any options for 2d or 1d periodic?
|
||||||
|
# allow for box size changes, e.g. every step for NPT
|
||||||
|
# check NWC func call error returns ?
|
||||||
|
|
||||||
|
import sys,time
|
||||||
|
from ctypes import *
|
||||||
|
|
||||||
|
import numpy as np
|
||||||
|
from mpi4py import MPI
|
||||||
|
import MDI_Library as mdi
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
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',
|
||||||
|
]
|
||||||
|
|
||||||
|
# atomic_number_to_symbol converts atomic number to element symbol
|
||||||
|
|
||||||
|
atomic_number_to_symbol = {}
|
||||||
|
for i,symbol in enumerate(ELEMENTS):
|
||||||
|
atomic_number_to_symbol[i+1] = symbol
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# global data
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
world = 0
|
||||||
|
me = nprocs = 0
|
||||||
|
|
||||||
|
if MPI._sizeof(MPI.Comm) == sizeof(c_int):
|
||||||
|
MPI_Comm = c_int
|
||||||
|
else:
|
||||||
|
MPI_Comm = c_void_p
|
||||||
|
|
||||||
|
exitflag = False
|
||||||
|
|
||||||
|
AIMD = 0
|
||||||
|
QMMM = 1
|
||||||
|
mode = AIMD
|
||||||
|
|
||||||
|
# NWChem PWDFT library
|
||||||
|
|
||||||
|
libname = "libpwdft.so"
|
||||||
|
libpwdft = None
|
||||||
|
|
||||||
|
# QM inputs
|
||||||
|
|
||||||
|
nw_template = ""
|
||||||
|
nw_infile = ""
|
||||||
|
nw_outfile = ""
|
||||||
|
|
||||||
|
flag_qm_natoms = flag_mm_natoms = 0
|
||||||
|
flag_box = flag_box_displs = 0
|
||||||
|
flag_qm_elements = 0
|
||||||
|
flag_qm_coords = flag_qm_potential = 0
|
||||||
|
flag_mm_elements = 0
|
||||||
|
flag_mm_coords = flag_mm_charges = 0
|
||||||
|
|
||||||
|
box = np.empty(9)
|
||||||
|
box_displ = np.empty(3)
|
||||||
|
|
||||||
|
qm_natoms = 0
|
||||||
|
qm_elements = None
|
||||||
|
qm_coords = None
|
||||||
|
qm_potential = None
|
||||||
|
|
||||||
|
mm_natoms = 0
|
||||||
|
mm_coords = None
|
||||||
|
mm_charges = None
|
||||||
|
mm_elements = None
|
||||||
|
|
||||||
|
# QM outputs
|
||||||
|
|
||||||
|
qm_pe = 0.0
|
||||||
|
qm_stress = np.empty(9)
|
||||||
|
qm_forces = None
|
||||||
|
qm_charges = None
|
||||||
|
|
||||||
|
mm_forces = None
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# print error message and halt
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def error(txt):
|
||||||
|
if me == 0: print("ERROR:",txt)
|
||||||
|
world.Abort()
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# process non-MDI options to PWDFT
|
||||||
|
# if this script is executed independently:
|
||||||
|
# args = command-line args
|
||||||
|
# if this script is invoked as a plugin library:
|
||||||
|
# args = passed via MDI
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def options(other_options):
|
||||||
|
global nw_template,nw_infile,nw_outfile
|
||||||
|
if len(other_options) != 3:
|
||||||
|
error("Invalid args to NWChem wrapper: template_file infile outfile")
|
||||||
|
nw_template = other_options[0]
|
||||||
|
nw_infile = other_options[1]
|
||||||
|
nw_outfile = other_options[2]
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# operate as an engine
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def mdi_engine(other_options):
|
||||||
|
global world,MPI_Comm,libpwdft
|
||||||
|
|
||||||
|
# get the MPI intra-communicator for this engine
|
||||||
|
|
||||||
|
world = mdi.MDI_MPI_get_world_comm()
|
||||||
|
me = world.Get_rank()
|
||||||
|
nprocs = world.Get_size()
|
||||||
|
|
||||||
|
# process non-MDI command line args
|
||||||
|
|
||||||
|
options(other_options)
|
||||||
|
|
||||||
|
# confirm PWDFT is being run as an engine
|
||||||
|
|
||||||
|
role = mdi.MDI_Get_Role()
|
||||||
|
if not role == mdi.MDI_ENGINE:
|
||||||
|
error("Must run NWChem as an MDI engine")
|
||||||
|
|
||||||
|
# supported MDI commands
|
||||||
|
|
||||||
|
mdi.MDI_Register_Node("@DEFAULT")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT","EXIT")
|
||||||
|
|
||||||
|
# driver --> engine
|
||||||
|
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">NATOMS")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">CELL")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">CELL_DISPLS")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">ELEMENTS")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">COORDS")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">POTENTIAL_AT_NUCLEI")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">NLATTICE")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">LATTICE_ELEMENTS")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">CLATTICE")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">LATTICE")
|
||||||
|
|
||||||
|
# engine --> driver
|
||||||
|
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT","<PE")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT","<FORCES")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT",">LATTICE_FORCES")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT","<STRESS")
|
||||||
|
mdi.MDI_Register_Command("@DEFAULT","<CHARGES")
|
||||||
|
|
||||||
|
# load PWDFT lib and set ctypes signatures for function calls
|
||||||
|
|
||||||
|
pwdft_load()
|
||||||
|
|
||||||
|
# one-time operation to establish a connection with the driver
|
||||||
|
|
||||||
|
mdicomm = mdi.MDI_Accept_Communicator()
|
||||||
|
|
||||||
|
# set callback to execute_command
|
||||||
|
|
||||||
|
mdi.MDI_Set_Execute_Command_Func(execute_command,None)
|
||||||
|
|
||||||
|
# infinite loop to exchange messages with driver
|
||||||
|
# until EXIT command received
|
||||||
|
|
||||||
|
while not exitflag:
|
||||||
|
command = mdi.MDI_Recv_Command(mdicomm)
|
||||||
|
command = world.bcast(command,root=0)
|
||||||
|
execute_command(command,mdicomm,None)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# called in loop by mdi_engine()
|
||||||
|
# called internally from MDI_Recv_command() until EXIT received
|
||||||
|
# command = name of MDI command
|
||||||
|
# mdicomm = MDI communicator for all MDI commands
|
||||||
|
# object_ptr = ptr to data if necessary
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def execute_command(command,mdicomm,object_ptr):
|
||||||
|
global exitflag
|
||||||
|
|
||||||
|
# driver says done
|
||||||
|
|
||||||
|
if command == "EXIT":
|
||||||
|
exitflag = True
|
||||||
|
|
||||||
|
# MDI commands which setup quantum calculation
|
||||||
|
|
||||||
|
elif command == ">NATOMS":
|
||||||
|
receive_qm_natoms(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">CELL":
|
||||||
|
receive_box(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">CELL_DISPL":
|
||||||
|
receive_box_displ(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">ELEMENTS":
|
||||||
|
receive_qm_elements(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">COORDS":
|
||||||
|
receive_qm_coords(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">POTENTIAL_AT_NUCLEI":
|
||||||
|
receive_qm_potential(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">NLATTICE":
|
||||||
|
receive_mm_natoms(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">LATTICE_ELEMENTS":
|
||||||
|
receive_mm_elements(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">CLATTICE":
|
||||||
|
receive_mm_coords(mdicomm)
|
||||||
|
|
||||||
|
elif command == ">LATTICE":
|
||||||
|
receive_mm_charges(mdicomm)
|
||||||
|
|
||||||
|
# MDI commands which retreive quantum results
|
||||||
|
# each may also trigger the quantum calculation
|
||||||
|
|
||||||
|
elif command == "<PE":
|
||||||
|
evaluate()
|
||||||
|
ierr = mdi.MDI_Send(qm_pe,1,mdi.MDI_DOUBLE,mdicomm)
|
||||||
|
if ierr: error("MDI: <PE data")
|
||||||
|
|
||||||
|
elif command == "<FORCES":
|
||||||
|
evaluate()
|
||||||
|
ierr = mdi.MDI_Send(qm_forces,3*qm_natoms,mdi.MDI_DOUBLE,mdicomm)
|
||||||
|
if ierr: error("MDI: <FORCES data")
|
||||||
|
|
||||||
|
elif command == "<LATTICE_FORCES":
|
||||||
|
evaluate()
|
||||||
|
ierr = mdi.MDI_Send(mm_forces,3*mm_natoms,mdi.MDI_DOUBLE,mdicomm)
|
||||||
|
if ierr: error("MDI: <LATTICE_FORCES data")
|
||||||
|
|
||||||
|
elif command == "<STRESS":
|
||||||
|
evaluate()
|
||||||
|
ierr = mdi.MDI_Send(qm_stress,1,mdi.MDI_DOUBLE,mdicomm)
|
||||||
|
if ierr: error("MDI: <STRESS data")
|
||||||
|
|
||||||
|
elif command == "<CHARGES":
|
||||||
|
evaluate()
|
||||||
|
ierr = mdi.MDI_Send(qm_charges,qm_natoms,mdi.MDI_DOUBLE,mdicomm)
|
||||||
|
if ierr: error("MDI: <CHARGES data")
|
||||||
|
|
||||||
|
# unrecognized command
|
||||||
|
|
||||||
|
else: error("Unrecognized MDI command")
|
||||||
|
|
||||||
|
return 0
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive count of QM atoms from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_qm_natoms(mdicomm):
|
||||||
|
global flag_qm_natoms,qm_natoms
|
||||||
|
flag_qm_natoms = 1
|
||||||
|
|
||||||
|
qm_natoms = mdi.MDI_Recv(1,mdi.MDI_INT,mdicomm)
|
||||||
|
qm_natoms = world.bcast(qm_natoms,root=0)
|
||||||
|
allocate("qm")
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive 3 simulation box edge vectors from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_box(mdicomm):
|
||||||
|
global flag_box
|
||||||
|
flag_box = 1
|
||||||
|
|
||||||
|
ierr = mdi.MDI_Recv(9,mdi.MDI_DOUBLE,mdicomm,buf=box)
|
||||||
|
if ierr: error("MDI: >CELL data")
|
||||||
|
world.Bcast(box,root=0)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive simulation box displacement vector from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_box_displ(mdicomm):
|
||||||
|
global flag_box_displ
|
||||||
|
flag_box_displ = 1
|
||||||
|
|
||||||
|
ierr = mdi.MDI_Recv(3,mdi.MDI_DOUBLE,mdicomm,buf=box_displ)
|
||||||
|
if ierr: error("MDI: >CELL_DISPL data")
|
||||||
|
world.Bcast(box_displ,root=0)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive QM atom coords from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_qm_coords(mdicomm):
|
||||||
|
global flag_qm_coords
|
||||||
|
flag_qm_coords = 1
|
||||||
|
|
||||||
|
if not qm_natoms: error("Cannot MDI >COORDS if # of QM atoms = 0")
|
||||||
|
|
||||||
|
ierr = mdi.MDI_Recv(3*qm_natoms,mdi.MDI_DOUBLE,mdicomm,buf=qm_coords)
|
||||||
|
if ierr: error("MDI: >COORDS data")
|
||||||
|
world.Bcast(qm_coords,root=0)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive Coulomb potential at QM nuclei from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_qm_potential(mdicomm):
|
||||||
|
global flag_qm_potential
|
||||||
|
flag_qm_potential = 1
|
||||||
|
|
||||||
|
if not qm_natoms: error("Cannot MDI >POTENTIAL_AT_NUCLEI if # of QM atoms = 0")
|
||||||
|
|
||||||
|
ierr = mdi.MDI_Recv(qm_natoms,mdi.MDI_DOUBLE,mdicomm,buf=qm_potential)
|
||||||
|
if ierr: error("MDI: >POTENTIAL_AT_NUCLEI data")
|
||||||
|
world.Bcast(qm_potential,root=0)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive QM atomic numbers from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_qm_elements(mdicomm):
|
||||||
|
global flag_qm_elements
|
||||||
|
flag_qm_elements = 1
|
||||||
|
|
||||||
|
if not qm_natoms: error("Cannot MDI >ELEMENTS if # of QM atoms = 0")
|
||||||
|
|
||||||
|
ierr = mdi.MDI_Recv(qm_natoms,mdi.MDI_INT,mdicomm,buf=qm_elements)
|
||||||
|
if ierr: error("MDI: >ELEMENTS data")
|
||||||
|
world.Bcast(qm_elements,root=0)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive count of MM atoms from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_mm_natoms(mdicomm):
|
||||||
|
global flag_mm_natoms,mm_natoms
|
||||||
|
flag_mm_natoms = 1
|
||||||
|
|
||||||
|
mm_natoms = mdi.MDI_Recv(1,mdi.MDI_INT,mdicomm)
|
||||||
|
mm_natoms = world.bcast(mm_natoms,root=0)
|
||||||
|
allocate("mm")
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive MM atomic numbers from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_mm_elements(mdicomm):
|
||||||
|
global flag_mm_elements
|
||||||
|
flag_mm_elements = 1
|
||||||
|
|
||||||
|
if not mm_natoms: error("Cannot MDI >LATTICE_ELEMENTS if # of MM atoms = 0")
|
||||||
|
|
||||||
|
ierr = mdi.MDI_Recv(mm_natoms,mdi.MDI_INT,mdicomm,buf=mm_elements)
|
||||||
|
if ierr: error("MDI: >LATTICE_ELEMENTS data")
|
||||||
|
world.Bcast(mm_elements,root=0)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive MM atom coords from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_mm_coords(mdicomm):
|
||||||
|
global flag_mm_coords
|
||||||
|
flag_mm_coords = 1
|
||||||
|
|
||||||
|
if not mm_natoms: error("Cannot MDI >CLATTICE if # of MM atoms = 0")
|
||||||
|
|
||||||
|
ierr = mdi.MDI_Recv(3*mm_natoms,mdi.MDI_DOUBLE,mdicomm,buf=mm_coords)
|
||||||
|
if ierr: error("MDI: >CLATTICE data")
|
||||||
|
world.Bcast(mm_coords,root=0)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# receive charge on MM atoms from driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def receive_mm_charges(mdicomm):
|
||||||
|
global flag_mm_charges
|
||||||
|
flag_mm_charges = 1
|
||||||
|
|
||||||
|
if not mm_natoms: error("Cannot MDI >LATTICE if # of MM atoms = 0")
|
||||||
|
|
||||||
|
ierr = mdi.MDI_Recv(mm_natoms,mdi.MDI_DOUBLE,mdicomm,buf=mm_charges)
|
||||||
|
if ierr: error("MDI: >LATTICE data")
|
||||||
|
world.Bcast(mm_charges,root=0)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# allocate persistent data for QM or MM atoms
|
||||||
|
# called when qm_natoms or mm_natoms is reset by MDI driver
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def allocate(which):
|
||||||
|
global qm_elements,qm_coords,qm_potential,qm_forces,qm_charges
|
||||||
|
global mm_elements,mm_coords,mm_charges,mm_forces
|
||||||
|
|
||||||
|
if which == "qm":
|
||||||
|
n = qm_natoms
|
||||||
|
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)
|
||||||
|
|
||||||
|
if which == "mm":
|
||||||
|
n = mm_natoms
|
||||||
|
mm_elements = np.empty(n,dtype=np.int32)
|
||||||
|
mm_coords = np.empty((n,3))
|
||||||
|
mm_charges = np.empty(n)
|
||||||
|
mm_forces = np.empty((n,3))
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# perform a quantum calculation via NWChem PWDFT
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def evaluate():
|
||||||
|
global mode
|
||||||
|
global flag_qm_natoms,flag_mm_natoms
|
||||||
|
global flag_box,flag_box_displs
|
||||||
|
global flag_qm_elements,flag_qm_coords,flag_qm_potential
|
||||||
|
global flag_mm_elements,flag_mm_coords,flag_mm_charges
|
||||||
|
global qm_pe,qm_stress,qm_forces,qm_charges
|
||||||
|
global mm_forces
|
||||||
|
global dm_previous
|
||||||
|
|
||||||
|
# just return if the QM system was already evaluated
|
||||||
|
# happens when multiple results are requested by driver
|
||||||
|
|
||||||
|
any_flag = flag_qm_natoms + flag_mm_natoms + flag_box + flag_box_displs + \
|
||||||
|
flag_qm_elements + flag_qm_coords + flag_qm_potential + \
|
||||||
|
flag_mm_elements + flag_mm_coords + flag_mm_charges
|
||||||
|
if not any_flag: return
|
||||||
|
|
||||||
|
# if any of these MDI commands received from LAMMPS,
|
||||||
|
# treat it as a brand new system
|
||||||
|
|
||||||
|
new_system = 0
|
||||||
|
if flag_qm_natoms or flag_mm_natoms: new_system = 1
|
||||||
|
if flag_qm_elements or flag_mm_elements: new_system = 1
|
||||||
|
if new_system:
|
||||||
|
if flag_mm_natoms or flag_qm_potential: mode = QMMM
|
||||||
|
else: mode = AIMD
|
||||||
|
|
||||||
|
# if new system, error check that all needed MDI calls have been made
|
||||||
|
|
||||||
|
if new_system:
|
||||||
|
if not flag_qm_natoms: error("QM atom count not specified")
|
||||||
|
if not flag_qm_elements or not flag_qm_coords or not flag_qm_potential:
|
||||||
|
error("QM atom properties not fully specified")
|
||||||
|
|
||||||
|
# setup new system within PWDFT
|
||||||
|
|
||||||
|
if new_system: pwdft_initialize()
|
||||||
|
|
||||||
|
# QMMM with QM and MM atoms
|
||||||
|
|
||||||
|
world_ptr = MPI._addressof(world)
|
||||||
|
c_world = MPI_Comm.from_address(world_ptr)
|
||||||
|
c_qm_pe = c_double(qm_pe)
|
||||||
|
c_nw_outfile = nw_outfile.encode()
|
||||||
|
|
||||||
|
if mode == QMMM:
|
||||||
|
print("QMMM minimizer")
|
||||||
|
time1 = time.time()
|
||||||
|
nwerr = libpwdft.\
|
||||||
|
c_lammps_pspw_qmmm_minimizer_filename(c_world,qm_coords,qm_potential,
|
||||||
|
qm_forces,qm_charges,byref(c_qm_pe),
|
||||||
|
False,True,c_nw_outfile)
|
||||||
|
# NOTE: check nwerr return?
|
||||||
|
qm_pe = c_qm_pe.value
|
||||||
|
time2 = time.time()
|
||||||
|
print("DONE QMMM minimizer",nwerr,time2-time1)
|
||||||
|
print("PE",qm_pe)
|
||||||
|
print("FORCE",qm_forces)
|
||||||
|
print("CHARGES",qm_charges)
|
||||||
|
|
||||||
|
# AIMD with only QM atoms
|
||||||
|
|
||||||
|
elif mode == AIMD:
|
||||||
|
print("AIMD minimizer")
|
||||||
|
time1 = time.time()
|
||||||
|
nwerr = libpwdft.\
|
||||||
|
c_lammps_pspw_aimd_minimizer_filename(c_world,qm_coords,qm_forces,
|
||||||
|
byref(c_qm_pe),c_nw_outfile)
|
||||||
|
# NOTE: check nwerr return?
|
||||||
|
qm_pe = c_qm_pe.value
|
||||||
|
time2 = time.time()
|
||||||
|
print("DONE AIMD minimizer",nwerr,time2-time1)
|
||||||
|
|
||||||
|
# clear flags for all MDI commands for next QM evaluation
|
||||||
|
|
||||||
|
flag_qm_natoms = flag_mm_natoms = 0
|
||||||
|
flag_box = flag_box_displs = 0
|
||||||
|
flag_qm_elements = 0
|
||||||
|
flag_qm_coords = flag_qm_potential = 0
|
||||||
|
flag_mm_elements = 0
|
||||||
|
flag_mm_coords = flag_mm_charges = 0
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# load PWDFT lib and set ctypes signatures for function calls
|
||||||
|
# set ctypes signatures for 3 function calls to PWDFT lib
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def pwdft_load():
|
||||||
|
global libpwdft
|
||||||
|
|
||||||
|
libpwdft = CDLL(libname,RTLD_GLOBAL)
|
||||||
|
|
||||||
|
libpwdft.c_lammps_pspw_input_filename.restype = None
|
||||||
|
libpwdft.c_lammps_pspw_input_filename.argtypes = \
|
||||||
|
[MPI_Comm, c_char_p, c_char_p]
|
||||||
|
|
||||||
|
nparray = np.ctypeslib.ndpointer(dtype=np.float64,ndim=2,flags="C_CONTIGUOUS")
|
||||||
|
npvector = np.ctypeslib.ndpointer(dtype=np.float64,ndim=1,flags="C_CONTIGUOUS")
|
||||||
|
|
||||||
|
libpwdft.c_lammps_pspw_qmmm_minimizer_filename.restype = c_int
|
||||||
|
libpwdft.c_lammps_pspw_qmmm_minimizer_filename.argtypes = \
|
||||||
|
[MPI_Comm, nparray, npvector, nparray, npvector, POINTER(c_double),
|
||||||
|
c_bool, c_bool, c_char_p]
|
||||||
|
|
||||||
|
libpwdft.c_lammps_pspw_aimd_minimizer_filename.restype = c_int
|
||||||
|
libpwdft.c_lammps_pspw_aimd_minimizer_filename.argtypes = \
|
||||||
|
[MPI_Comm, nparray, nparray, POINTER(c_double), c_char_p]
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# create PWDFT input file with box and list of atoms
|
||||||
|
# invoke PWDFT function to read it
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def pwdft_initialize():
|
||||||
|
|
||||||
|
# box, qm_coords, mm_coords must be converted to Angstroms
|
||||||
|
|
||||||
|
angstrom_to_bohr = mdi.MDI_Conversion_factor("angstrom","bohr")
|
||||||
|
bohr_to_angstrom = 1.0 / angstrom_to_bohr
|
||||||
|
|
||||||
|
box_A = box * bohr_to_angstrom
|
||||||
|
qm_coords_A = qm_coords * bohr_to_angstrom
|
||||||
|
|
||||||
|
# proc 0 reads template file, writes PWDFT input file
|
||||||
|
|
||||||
|
if me == 0:
|
||||||
|
lines = open(nw_template,'r').readlines()
|
||||||
|
|
||||||
|
fp = open(nw_infile,'w')
|
||||||
|
|
||||||
|
for line in lines:
|
||||||
|
word = line.strip()
|
||||||
|
if word == "GEOMINSERT":
|
||||||
|
print("geometry noautosym noautoz nocenter",file=fp);
|
||||||
|
print("system crystal cartesian",file=fp)
|
||||||
|
print("lattice_vectors",file=fp)
|
||||||
|
print("%20.16g %20.16g %20.16g" % (box_A[0],box_A[1],box_A[2]),file=fp)
|
||||||
|
print("%20.16g %20.16g %20.16g" % (box_A[3],box_A[4],box_A[5]),file=fp)
|
||||||
|
print("%20.16g %20.16g %20.16g" % (box_A[6],box_A[7],box_A[8]),file=fp)
|
||||||
|
print("end\n",file=fp)
|
||||||
|
|
||||||
|
for i in range(qm_natoms):
|
||||||
|
symbol = atomic_number_to_symbol[qm_elements[i]]
|
||||||
|
print("%s %20.16g %20.16g %20.16g" %
|
||||||
|
(symbol,qm_coords_A[i][0],qm_coords_A[i][1],qm_coords_A[i][2]),
|
||||||
|
file=fp)
|
||||||
|
print("end\n",file=fp)
|
||||||
|
|
||||||
|
else: print(line,file=fp,end="")
|
||||||
|
|
||||||
|
fp.close()
|
||||||
|
|
||||||
|
# all procs call pspw_input_filename() which processes input file
|
||||||
|
# performs initial QM calculation within PWDFT
|
||||||
|
|
||||||
|
world_ptr = MPI._addressof(world)
|
||||||
|
c_world = MPI_Comm.from_address(world_ptr)
|
||||||
|
infile = nw_infile.encode()
|
||||||
|
outfile = nw_outfile.encode()
|
||||||
|
|
||||||
|
print("INPUT filename")
|
||||||
|
time1 = time.time()
|
||||||
|
nwerr = libpwdft.c_lammps_pspw_input_filename(c_world,infile,outfile)
|
||||||
|
time2 = time.time()
|
||||||
|
print("DONE INPUT filename",nwerr,time2-time1)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# function called by MDI driver
|
||||||
|
# only when it invokes pyscf_mdi.py as a plugin
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
def MDI_Plugin_init_nwchem_mdi(plugin_state):
|
||||||
|
|
||||||
|
# other_options = all non-MDI args
|
||||||
|
# -mdi arg is processed and stripped internally by MDI
|
||||||
|
|
||||||
|
other_options = []
|
||||||
|
|
||||||
|
mdi.MDI_Set_plugin_state(plugin_state)
|
||||||
|
narg = mdi.MDI_Plugin_get_argc()
|
||||||
|
|
||||||
|
for iarg in range(narg):
|
||||||
|
arg = mdi.MDI_Plugin_get_arg(iarg)
|
||||||
|
other_options.append(arg)
|
||||||
|
|
||||||
|
# start running as an MDI engine
|
||||||
|
|
||||||
|
mdi_engine(other_options)
|
||||||
|
|
||||||
|
# --------------------------------------------
|
||||||
|
# main program
|
||||||
|
# invoked when MDI driver and pyscf_mdi.py
|
||||||
|
# are run as independent executables
|
||||||
|
# --------------------------------------------
|
||||||
|
|
||||||
|
if __name__== "__main__":
|
||||||
|
|
||||||
|
# mdi_option = single arg in quotes that follows -mdi
|
||||||
|
# other_options = all non-MDI args
|
||||||
|
|
||||||
|
mdi_option = ""
|
||||||
|
other_options = []
|
||||||
|
|
||||||
|
narg = len(sys.argv)
|
||||||
|
args = sys.argv
|
||||||
|
|
||||||
|
iarg = 1
|
||||||
|
while iarg < narg:
|
||||||
|
arg = args[iarg]
|
||||||
|
if arg == "-mdi" or arg == "--mdi":
|
||||||
|
if narg > iarg+1: mdi_option = sys.argv[iarg+1]
|
||||||
|
else: error("NWChem -mdi argument not provided")
|
||||||
|
iarg += 1
|
||||||
|
else: other_options.append(arg)
|
||||||
|
iarg += 1
|
||||||
|
|
||||||
|
if not mdi_option: error("NWChem -mdi option not provided")
|
||||||
|
|
||||||
|
# call MDI_Init with just -mdi option
|
||||||
|
|
||||||
|
mdi.MDI_Init(mdi_option)
|
||||||
|
|
||||||
|
# start running as an MDI engine
|
||||||
|
|
||||||
|
mdi_engine(other_options)
|
||||||
17
examples/QUANTUM/NWChem/template.water.nw
Normal file
17
examples/QUANTUM/NWChem/template.water.nw
Normal file
@ -0,0 +1,17 @@
|
|||||||
|
Title "LAMMPS wrapping of PWDFT"
|
||||||
|
|
||||||
|
memory 1900 mb
|
||||||
|
|
||||||
|
echo
|
||||||
|
|
||||||
|
GEOMINSERT
|
||||||
|
|
||||||
|
nwpw
|
||||||
|
xc pbe
|
||||||
|
cutoff 30.0
|
||||||
|
2d-hcurve
|
||||||
|
tolerances 1.0e-9 1.0-9
|
||||||
|
apc on
|
||||||
|
end
|
||||||
|
|
||||||
|
task pspw gradient
|
||||||
@ -1,5 +1,5 @@
|
|||||||
Each of the directories shows how to use LAMMPS in tandem with a
|
Each directory shows how to use LAMMPS in tandem with a specific
|
||||||
specific quantum code
|
quantum code:
|
||||||
|
|
||||||
LATTE = semi-empirical tight-binding code from LANL
|
LATTE = semi-empirical tight-binding code from LANL
|
||||||
https://www.osti.gov/biblio/1526907-los-alamos-transferable-tight-binding-energetics-latte-version
|
https://www.osti.gov/biblio/1526907-los-alamos-transferable-tight-binding-energetics-latte-version
|
||||||
@ -8,14 +8,14 @@ LATTE = semi-empirical tight-binding code from LANL
|
|||||||
PySCF = ??? from Caltech
|
PySCF = ??? from Caltech
|
||||||
add link
|
add link
|
||||||
|
|
||||||
|
NWChem = computational chemistry code from PNNL
|
||||||
|
focus here is on DFT portion of NWChem = PWDFT library
|
||||||
|
https://www.nwchem-sw.org
|
||||||
|
|
||||||
-----------------------------------------------------
|
-----------------------------------------------------
|
||||||
|
|
||||||
To be added later (as of Jan 2023):
|
To be added later (as of Jan 2023):
|
||||||
|
|
||||||
NWChem = computational chemistry code from PNNL
|
|
||||||
focus here is on DFT portion of NWChem = PWDFT
|
|
||||||
https://www.nwchem-sw.org
|
|
||||||
|
|
||||||
Quantum Espresso (QE) = DFT code for materials modeling
|
Quantum Espresso (QE) = DFT code for materials modeling
|
||||||
https://www.quantum-espresso.org/
|
https://www.quantum-espresso.org/
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user